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Introduction 


An  important  initial  screening  step  in  the  detection  of  breast  cancer  is  the  ability  to  identify  select 
areas  of  atypical  density  that  require  further  evaluation.  Currently,  mammography  is  the  clinical 
standard  for  screening  and  provides  useful  but  at  times  ambiguous  information,  which  can 
necessitate  further  invasive  workup  of  benign  lesions.  Alternative  methods  such  as  elastography 
have  shown  potential  in  enhancing  the  diagnostic  process  by  providing  information  about  the 
tissue  composition  [1,2].  Modality-independent  elastography  (MIE)  is  a  novel  image  processing 
technique  that  combines  finite  element  models  of  soft-tissue  deformation  with  measures  of  image 
similarity  in  order  to  reconstruct  elastic  property  distributions  throughout  the  tissue.  The  basic 
requirements  for  the  method  are  two  images  of  the  tissue  in  different  states  of  deformation  (e.g. 
compression).  MIE  then  updates  the  estimate  of  the  material  properties  via  a  matching  process 
between  the  two  images.  The  final  result  is  a  map  of  the  breast  (or  other  tissue  of  interest)  that 
reflects  material  inhomogeneity,  such  as  in  the  case  of  a  tumor  mass  that  disrupts  the  surrounding 
structure  of  normal  tissue.  Because  MIE  works  on  probing  the  differences  between  images,  it 
can  be  used  to  not  only  work  in  concert  with  more  traditional  screening  techniques  but  also 
address  a  possible  gap  when  those  methods  are  unable  to  directly  discern  tissues  of  interest. 


Body 

As  stated  in  the  original  proposal,  three  main  aims  of  this  project  are  to  (1)  expand  and  refine  the 
current  MIE  technique  to  enhance  its  efficiency  and  capabilities,  (2)  to  perform  analyses  on 
texture  in  input  images  and  quantify  statistical  parameters  capable  of  estimating  and  evaluating 
the  success  of  elastographic  reconstruction,  and  (3)  to  engineer  a  device  that  can  accurately 
produce  compressive  forces  necessary  for  phantom  setups  within  current  imaging  systems, 
providing  the  basis  for  a  future  device  that  can  be  used  in  a  clinical  setting.  In  this  past  year, 
progress  on  all  three  aims  has  been  made.  The  original  specific  aim  and  the  relevant  proposed 
work  for  each  is  listed  below  and  addressed. 

Specific  Aim  #1  stated:  “To  expand  and  refine  the  current  MIE  technique  to  enhance  its 
efficiency,  as  well  as  add  new  capabilities  such  as  handling  a  full  3D  or  combined  2D/3D 
elastodynamic  model  for  improved  accuracy.” 

A  framework  utilizing  parallel  processing  techniques  has  been  designed  and  implemented  to 
accommodate  a  fully  three-dimensional  version  of  the  MIE  algorithm.  A  key  component,  but 
also  the  bottleneck,  of  the  MIE  methodology  is  the  calculation  of  a  Jacobian  matrix  sensitive  to 
spatial  discretization.  Each  column  of  the  Jacobian  matrix  reprsents  a  single  “forward  solve”  of 
the  problem  and  corresponds  to  a  single  ‘elasticity  region’  within  the  domain,  encompassing  a 
solution  of  the  finite  element  model,  image  deformation,  and  similarity  metric  calculation.  In 
2D,  a  fairly  large  number  of  regions  (about  400)  were  already  needed  to  adequately  sample  the 
area  of  study  and  attempt  to  search  for  a  presumed  lesion.  For  an  organ  that  is  most  properly 
analyzed  as  a  three-dimensional  system,  such  as  the  breast,  our  initial  empire  calculations 
indicated  that  the  increase  in  dimensionality  would  necessitate  upwards  of  3200  regions  to  be 
partitioned  over  a  [typical]  volume.  In  this  case,  the  size  of  the  finite  element  model  itself  is 
increased  by  nearly  two  orders  of  magnitude  as  well.  It  was  therefore  apparent  that  the  code 
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complexity  and  computational  resource  demands  of  3D  MIE  far  surpass  the  capabilities  of  the 
original  2D  MATLAB/FORTRAN/LAPACK  design.  In  order  to  perform  what  would  otherwise 
be  an  intractable  computational  problem,  the  MPI  standard  of  parallel  communication  was  used 
to  split  up  the  Jacobian  formation  process  as  a  static  SPMD  (single-process,  multiple-data) 
scheme  within  a  C/C++  Gauss-Newton  optimization.  In  addition,  the  Portable  Extensible 
Toolkit  for  Scientific  Computation  (PETSc)  toolkit  [3,4]  was  interfaced  to  provide  the  necessary 
efficient  sparse  matrix  system  solvers.  We  utilized  a  share  of  100  CPUs  on  the  Vanderbilt 
University  Advanced  Computing  Center  for  Research  and  Education  cluster  to  demonstrate  that 
each  iteration  of  the  optimization  could  be  achieved  on  the  order  of  30  minutes,  as  opposed  to 
original  estimates  on  available  sequential  processing  machines  of  upwards  of  5000  minutes.  We 
note,  however,  that  to  effectively  traverse  the  entire  multi-dimensional  objective  function  space 
requires  several  (perhaps  tens  of)  iterations  would  be  required,  underscoring  the  innately  high 
computational  load  of  the  method. 

We  have  completed  a  preliminary  study  that  demonstrates  the  new  MIE  framework  in  action, 
using  simulated  deformation  based  on  clinical  data  obtained  from  both  CT  and  MR  scans.  The 
test  cases  involved  the  simulated  implantation  of  a  stiff  and  radiographically  occult  2-cm 
spherical  tumor  at  the  center  of  the  breast.  Guided  by  a  finite  element  mesh  deformation  using 
prescribed  boundary  displacements  (designed  to  mimic  a  compression  source  as  described  in 
Specific  Aim  3),  a  target  image  volume  was  created.  Discretizations  of  the  model  and  image 
domains  as  per  the  MIE  methodology  were  then  used  to  reconstruct  the  inclusion.  In  both  cases, 
MIE  indicated  that  a  lesion  ~2x  stiffer  than  surrounding  tissue  existed  and  was  localized  to  the 
known  position.  However,  due  to  inexact  partitioning  of  the  mesh  elements,  the  algorithm 
actually  detected  a  tumor  closer  to  3  cm  in  diameter,  leading  to  a  pseudo-compensatory  decrease 
in  elasticity  contrast  in  the  model  (the  true  difference  in  stiffness  was  a  factor  of  6).  This  was 
explained  by  performing  an  a  priori  classification  of  the  mesh  elements  into  two  material  types 
according  to  their  spatial  membership  in  the  3200  regions.  After  this  adjustment  in  overlap,  the 
constrained  reconstructions  much  more  favorably  followed  an  objective  function  minimization 
that  resulted  in  a  stiffness  contrast  matching  the  original  reconstruction.  Based  on  these 
experiments,  it  is  our  current  assertion  that  shaping  the  objective  function  by  dynamically 
rearranging  the  spatial  discretization  of  the  model  during  the  optimization  can  lead  to  improved 
elasticity  contrast  resolution,  and  studies  are  underway  to  address  this  issue.  The  details  of  the 
work  described  are  provided  within  Appendix  A  of  this  document. 

As  a  further  note,  we  acknowledge  that  originially  proposed  work  involving  the  use  of  a 
combined  2D-3D  model  is  still  being  considered  as  a  means  of  reducing  the  scope  of  the  problem 
with  respect  to  dimsenionality.  However,  the  current  implementation  is  deemed  viable  with 
sufficient  allocation  of  computational  resources.  In  addition,  further  research  into  the  effects  of 
dimensionality  was  pursued  in  the  context  of  another  possible  application  for  MIE,  that  is,  in 
dermoscopic  examination.  We  simulated  a  spectrum  of  geometrical  variation  in  size  for  a 
posited  lesion  embedded  within  skin  (e.g.  radius  and  depth  of  penetration).  Our  findings  indicate 
that  while  a  thin-membrane  system  is  capable  of  differentiating  a  stiff  object  from  its 
surroundings  for  most  scenarios,  inevitably,  a  2D  approximation  of  the  model  has  its  limitations, 
and  it  is  quite  difficult  to  avoid  the  need  to  account  for  full  spatial  dimensionality.  For  further 
discussion,  please  refer  to  Appendix  B  of  this  document. 
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Specific  Aim  #2  stated:  “To  perform  texture  analysis  on  input  images  in  order  to  quantify  a 
statistical  parameter  capable  of  estimating  the  success  of  elastographic  reconstruction.” 

Tolerance  to  improper  input  has  been  tested  with  statistical  quantification  of  reconstruction 
success.  Our  observations  during  the  ongoing  development  and  testing  of  the  MIE  method 
prompted  questions  concerning  the  quality  of  data  necessary  and  sufficient  to  achieve 
satisfactory  results,  as  well  as  how  to  evaluate  the  fidelity  of  the  reconstructed  elasticity  image. 
The  primary  inputs  to  the  reconstruction  method  are  the  acquired  images  and  the  delineated 
boundary  conditions  on  the  domain.  Our  investigation  began  by  examining  a  two-dimensional 
system  and  attempting  to  defeat  the  algorithm  with  increasing  levels  of  degraded  input.  We 
briefly  summarize  the  key  results  of  these  studies  as  described  in  a  previous  Annual  Report.  The 
first  experiment  used  additive  image  noise  to  obscure  the  underlying  texture  to  reflect  possible 
scenarios  of  corruption  during  acquisition.  Noise  fields  were  created  from  a  zero-mean  Gaussian 
random  distribution  along  the  variance  of  non-background  pixels  and  scaled  according  to  the 
total  power  at  1,  5,  10,  15,  20,  25,  and  30%.  It  was  found  that  the  reconstruction  was  tolerant  of 
image  noise  up  to  approximately  10%.  Figure  1  demonstrates  the  degradative  effects  of  image 
noise.  The  second  experiment  involved  boundary  condition  selection  error.  A  gold  standard  set 
of  boundary  conditions  known  to  produce  an  accurate  reconstruction  was  modified  using 
randomized  vectors  of  equal  magnitude  (0.1,  0.2,  0.3,  0.5,  0.75,  1.0,  1.5,  and  2.0  pixel  units) 
reflecting  a  range  of  typical  localization  skill  for  users  from  poor  to  expert.  It  was  noted  that 
randomizing  all  the  vectors  can  actually  result  in  twisting  of  elements  that  resulted  in  significant 
alterations  of  displacements  in  the  interior  of  the  mesh,  leading  to  grossly  inaccurate  model 
deformations. 


Figure  1.  2D  reconstructions  resulting  from  the  distortion  of  the  target  image  using  additive  Gaussian  random  noise 
(from  top  left:  1,  5,  10,  20,  25,  30%).  The  true  elasticity  distribution  is  a  centrally  located  and  roughly  circular 
region,  and  the  noise  progressively  confounds  the  reconstruction. 


We  also  note  that  the  ‘quality  of  reconstruction  score’  (QRS)  was  developed  within  the  context 
of  this  project  in  order  to  quantify  the  retrospective  localization  accuracy  of  the  method.  In 
comparison  and  conjuction  with  more  standard  measures  in  the  elastography  field  such  as 
contrast-to-noise  ratio  [5],  QRS  has  been  used  to  determine  relevant  positional  and  material 
characterization  in  both  simulation  and  data  studes.  The  metric  is  determined  by  a  classification 
of  the  reconstruction  [6]  that  is  then  compared  to  the  (known)  segmentation  of  the  actual 


3 


elasticity  distribution.  By  examining  the  rate  of  accurately  selecting  an  element  of  the  domain  to 
be  of  the  correct  material  class,  a  conditional  probability  closely  related  to  the  positive  predictive 
value  of  the  test  is  obtained;  we  have  determined  a  posteriori  that  a  QRS>80%  is  typically 
indicative  of  a  successful  reconstruction.  The  use  of  QRS  was  similarly  applied  to  the  analysis  of 
our  3D  MIE  simulation  reconstuction  experiments  as  included  in  Appendix  A. 

While  we  are  continuing  to  examine  the  possibility  of  analyzing  image  texture  prospectively  for 
the  purposes  of  predicting  MIE  success,  we  acknowledge  that  variability  among  acquisition 
systems  and  subjects  has  thus  far  made  this  extremely  difficult.  Proposed  work  involving  the  use 
of  a  feature  tracking  and  frequency  domain  analysis  is  under  investigation  but  not  completed  at 
this  time.  As  more  data  is  collected,  it  is  hoped  that  establishing  a  pattern  for  understanding  the 
reconstruction  algorithm  behavior  will  become  statistically  relevant.  However,  partially  based 
on  our  preliminary  findings,  we  believe  that  the  natural  heterogeneity  of  tissue  is  typically 
sufficient  for  reconstruction  purposes.  Our  efforts  have  therefore  been  refocused  on  addressing 
the  practical  aspects  of  implementing  MIE  outside  of  an  in  silico  environment,  especially  with 
regards  to  the  degradation  of  quality  due  to  the  (mis)estimation  in  the  process  of  boundary 
condition  assignment.  It  is  the  propagation  of  error  from  this  source  into  the  model  and 
downstream  to  the  other  components  of  the  algorithm  that  appears  to  be  the  greatest  contributor 
in  an  unsatisfactory  reconstruction  (as  demonstrated  in  the  work  presented  in  Appendix  A  and 
C).  In  pursuing  this  line  of  work,  a  number  of  novel  as  well  as  standard  means  for  determining 
point  correspondence  to  provide  an  accurate  surface  registration  have  been  developed  and  tested. 
At  this  time,  our  conclusion  is  that  utilizing  a  thin-plate  spline  is  the  best  method.  Because  this 
requires  the  tracking  of  surface  markers,  an  level  set-matching  method  using  Laplace’s  equation 
is  deemed  appropriate  for  use  in  the  event  that  such  fiducials  are  not  available. 

Specific  Aim  #3  stated:  “To  engineer  a  device  that  can  accurately  produce  compressive 
forces  necessary  for  phantom  setups  within  current  clinical  imaging  systems,  providing  the 
basis  for  a  future  device  that  can  be  used  in  a  clinical  setting.” 

A  compression  device  has  been  constructed  and  tested  in  magnetic  resonance  (MR)  and  X-ray 
computed  tomography  (CT)  imaging  systems  using  a  polyvinyl  alcohol  phantom  and  contrast 
agents.  The  compression  device  is  composed  of  a  rectangular  Plexiglas  frame  that  traps  the 
phantom  in  at  least  two  directions  with  a  sliding  wall  and  the  compression  plate,  which  houses  an 
air  bladder  in  a  polycarbonate  frame  (see  Figure  2).  When  inflated,  the  air  bladder  provides  a 
deformation  of  up  to  5  cm.  The  prototypical  phantom  used  has  been  fabricated  as  a  polyvinvyl 
alcohol  cryogel  (~250  cc,  7%  wt/vol)  in  a  manner  consistent  with  the  methods  presented  in  [7, 

8],  A  notable  difference  in  the  manufacturing  process  has  been  the  addition  of  10%  v/v  glycerol 
to  the  base  solution.  A  spherical  mold  is  used  to  create  a  tumor  replica  after  the  application  of  a 
24-hour  freeze-thaw  cycle,  which  is  then  suspended  within  a  primary  mold  for  the  phantom  and 
frozen  again  within  the  remaining  bulk  material.  Contrast  agents  (e.g.  iodine,  barium,  or  copper 
sulfate)  are  added  as  needed  to  provide  intensity  variation  depending  on  the  imaging  modality. 
The  result  is  a  domed  breast-like  shape  approximately  10-11  cm  at  the  base  and  tapering  over  a 
depth  of  about  5-6  cm.  A  study  of  this  system  as  used  in  a  CT  scanner  has  been  completed,  with 
elastic  contrasts  obtained  by  MIE  reconstruction  having  been  compared  against  independent 
mechanical  testing  of  the  cryogel  materials.  Key  results  obtained  are  being  prepared  for  a 
manuscript  to  be  submitted  soon  after  the  preparation  of  this  report. 
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Figure  2.  Top  row:  Photographs  of  the  polyvinyl  alcohol  phantom  inside  the  compression  device  without  (left)  and 
with  compression  (right). 


A  prototype  compression  chamber  that  is  more  clinically  oriented  has  been  designed  to  fit  into 
the  chassis  of  a  Philips  Intera  MR  breast  coil  unit.  It  has  been  fabricated  from  clear  acrylic  tube 
segments  in  which  the  air  bladders  are  attached  using  polycarbonate  pins  and  then  covered  with 
an  expandable  nylon  sheet.  A  small  series  of  human  volunteer  trials  have  been  performed  to  test 
the  usage  of  the  device.  To  this  date,  no  negative  feedback  has  been  given  regarding  its  effects 
on  the  patient.  Figure  3  below  shows  the  schematic  design  of  the  system  as  well  as  an  overhead 
view  of  its  assembly.  Further  image  acquisition  studies  are  currently  being  designed  to  utilize 
this  system. 


Figure  3.  Left:  Schematic  of  compression  device  designed  for  clinical  use  breast  coil.  Right:  Photograph  of 
assembly  looking  down  into  the  Philips  Intera  chassis. 
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Key  Research  Accomplishments 

•  3D  implementation  -  Computationally  complex  challenges  have  been  addressed  to  create  an 
effective  parallelized  version  of  MIE  for  full  three-dimensional  analysis 

•  Phantom  and  material  creation  -  design,  fabrication,  and  testing  of  PVA-C  for  the  creation  of 
appropriate  tissue  mimicking  systems 

•  Testing  with  phantoms  -  Initial  testing  has  been  performed  on  semi-anthropomorphic  and 
appropriate  breast  phantoms  for  elastography  research.  Initial  testing  has  also  been  performed  on 
phantoms  designed  for  surgical  simulations 


Reportable  Outcomes 

Work  on  the  MIE  method  during  this  period  has  resulted  in  three  conference  papers  and  an 
additional  poster  presentation  as  well  as  acceptance  of  a  peer-reviewed  journal  article.  Didactic 
coursework  requirements  for  the  PhD  degree  have  been  completed  at  this  time  as  well  as 
successful  presentation  of  the  thesis  proposal  to  a  faculty  committee. 

Poster  Presentations 
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Conclusions 


The  current  results  and  progress  denoted  in  this  report  are  within  the  proposed  statement  of  work 
and  are  encouraging  towards  completion  of  the  overall  objectives  with  further  effort.  No 
significant  deviations  are  reported  at  this  time,  although  in  certain  regards  as  noted  within  the 
text,  our  efforts  are  being  refocused  towards  relevant  practical  issues  of  the  project  that  have 
become  more  prominent  as  the  work  has  progressed.  Notable  accomplishments  at  this  time 
include  the  implementation  of  a  fully  3D  MIE  application,  further  analysis  of  factors  influencing 
elastographic  reconstruction  fidelity,  and  engineering  of  devices  suitable  for  MIE  data 
acquisition  in  both  phantom  and  human  studies.  Future  work  will  be  directed  towards  following 
up  in  all  of  these  areas,  as  well  as  newer  studies  into  ex  vivo  tissue  characterization  and 
collaborative  clinical  efforts  generated  from  interest  regarding  our  project. 
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Abstract 

This  paper  reports  on  the  development  and  preliminary  testing  of  a  three- 
dimensional  implementation  of  an  inverse  problem  technique  for  extracting 
soft-tissue  elasticity  information  via  non-rigid  model-based  image  registration. 
The  modality-independent  elastography  (MIE)  algorithm  adjusts  the  elastic 
properties  of  a  biomechanical  model  to  achieve  maximal  similarity  between 
images  acquired  under  different  states  of  static  loading.  A  series  of  simulation 
experiments  with  clinical  image  sets  of  human  breasts  were  performed  to  test 
the  ability  of  the  method  to  identify  and  characterize  a  radiographically  occult 
stiff  lesion.  Because  boundary  conditions  are  a  critical  input  to  the  algorithm,  a 
comparison  of  three  methods  for  semi-automated  surface  point  correspondence 
was  conducted  in  the  context  of  systematic  and  randomized  noise  processes. 
The  results  illustrate  that  3D  MIE  was  able  to  successfully  reconstruct  elasticity 
images  using  data  obtained  from  both  magnetic  resonance  and  x-ray  computed 
tomography  systems.  The  lesion  was  localized  correctly  in  ah  cases  and  its 
relative  elasticity  found  to  be  reasonably  close  to  the  true  values  (3.5%  with 
the  use  of  spatial  priors  and  11.6%  without).  In  addition,  the  inaccuracies 
of  surface  registration  performed  with  thin-plate  spline  interpolation  did  not 
exceed  empiric  thresholds  of  unacceptable  boundary  condition  error. 


1.  Introduction 

Breast  cancer  is  the  most  common  cancer  of  women  in  the  United  States,  the  second  most 
common  cause  of  cancer  death  in  women  and  the  leading  cause  of  death  in  women  aged  45 
to  55.  Estimates  for  the  year  2007  indicate  that  178  480  American  women  will  be  diagnosed 
with  the  disease  and  40  910  women  will  die  from  it  (ACS  2007).  While  many  advances  have 
been  made  in  the  treatment  of  the  disease,  the  ability  to  detect  its  presence  for  either  screening 
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or  diagnostic  purposes  remains  an  area  of  active  research  involving  many  novel  forms  of 
imaging.  The  characterization  of  the  mechanical  properties  of  breast  tissue  is  an  important 
potential  source  of  clinical  information  because  of  the  long-standing  association  of  palpable 
differences  in  stiffness  with  possible  pathological  states.  A  minimally  invasive  methodology 
for  analyzing  tissue  elasticity  through  imaging  and/or  image  processing  techniques  is  a  central 
goal  of  the  field  of  elastography  (Parker  et  al  2005),  with  the  application  of  various  techniques 
being  found  not  only  in  the  interrogation  of  the  breast  (McKnight  et  al  2002,  Melodelima 
et  al  2006,  Sinkus  et  al  2000),  but  also  in  skin  (Miga  et  al  2005,  Tsap  et  al  1998,  Zhang  et  al 
2004),  prostate  (Curiel  et  al  2005,  Egorov  et  al  2006)  and  other  accessible  organ  systems. 

Many  of  the  current  elastography  methods  are  founded  in  ultrasound  (US)  (Ophir  et  al 
1991,  2000)  and  magnetic  resonance  (MR)  (Manduca  et  al  2001,  Muthupillai  et  al  1995) 
imaging  and  involve  the  estimation  of  induced  displacements  within  the  tissue  of  interest  to 
infer  the  elasticity  distribution.  We  have  recast  the  problem  as  a  physically  constrained,  non- 
rigid  image  registration  utilizing  numerical  models  of  static  deformation  with  image  similarity 
metrics  to  reconstruct  the  spatial  distribution  of  elasticity  parameters.  This  technique  has 
been  termed  ‘modality-independent  elastography’  (MIE)  (Miga  2002,  2003,  Washington  and 
Miga  2004)  because  of  its  ability  to  handle  anatomical  images  from  different  sources  with 
relatively  simple  modifications  to  the  acquisition  procedure.  To  date,  data  from  MR,  x-ray 
computed  tomography  (CT)  and  digital  photography  have  been  used  to  successfully  drive  the 
algorithm  in  two-dimensional  (2D)  work.  Others  have  also  pursued  similar  approaches  within 
the  context  of  ultrasound  elastography  (Garra  et  al  1997,  Gokhale  et  al  2004,  Sarvazyan  et  al 
1995),  optical  image  analysis  (Tsap  et  al  1998)  and  to  a  lesser  extent  with  magnetic  resonance 
elastography  (Fowlkes  et  al  1995).  While  the  use  of  MIE  in  2D  has  been  illuminating  for 
algorithmic  development  and  may  have  its  own  applications  in  studying  the  more  planar 
system  of  the  skin,  ultimately,  translation  of  the  method  to  utilize  volumetric  data  is  desirable 
(if  not  necessary)  in  order  to  provide  an  accurate  representation  of  organs  such  as  the  breast  as  a 
whole.  In  this  work,  we  present  a  newly  realized  three-dimensional  (3D)  version  of  MIE  along 
with  simulation  experiments  to  evaluate  its  performance.  In  addition,  some  potential  effects  of 
degraded  input  quality  are  addressed  by  examining  robustness  of  the  algorithm  to  inaccuracies 
under  specified  boundary  conditions  and  then  comparing  the  reconstruction  fidelity  of  three 
different  techniques  developed  for  semi-automatic  generation  of  boundary  conditions. 

2.  Methods  and  materials 

2.7.  MIE  reconstruction  framework 

The  conceptual  framework  for  our  elastographic  reconstruction  has  been  previously  described 
in  Miga  (2002,  2003),  Miga  et  al  (2005)  and  Washington  and  Miga  (2004)).  To  review,  an 
image  of  a  tissue  of  interest  {source)  is  deformed  by  a  biomechanical  model  and  compared 
against  an  acquired  image  of  the  same  tissue  in  a  mechanically  loaded  state  {target).  Iterative 
updates  of  elasticity  parameters  to  the  model  are  performed  until  a  suitable  match  in  image 
similarity  is  achieved  in  a  least- squares  manner  to  satisfy  a  nonlinear  optimization  scheme. 
This  process  as  illustrated  in  figure  1  can  be  classified  as  an  inverse  problem,  with  model-based 
deformation  and  registration  of  the  source  image  representing  the  forward  problem. 

The  three  major  components  of  the  reconstruction  algorithm  are  the  biomechanical  model, 
image  comparison  and  the  optimization  routine.  Although  there  are  a  number  of  models  for 
soft- tissue  mechanics,  it  is  reasonably  appropriate  to  begin  with  a  general  elastic  body.  The 
partial  differential  equation  (PDEs)  that  expresses  a  state  of  mechanical  equilibrium  is 

V  •  cr  =  0  (1) 
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Figure  1.  Schematic  of  MIE  framework.  After  acquisition  of  image  data,  surface  representations 
are  segmented  from  the  pre-  and  post-deformation  volumes  ( source  and  target,  respectively). 
A  number  of  pre-processing  steps  are  performed  to  generate  boundary  conditions  for  the 
biomechanical  model,  which  produces  a  deformed  image  that  can  be  compared  with  the  true 
target  volume.  The  optimization  routine  updates  the  elasticity  distribution  until  the  best  similarity 
is  achieved. 


where  o  is  the  Cartesian  stress  tensor  (Boresi  and  Chong  1999).  We  have  elected  to  describe 
the  constitutive  tissue  behavior  using  Hooke’s  law  of  linear  elasticity,  which  states  that  the 
strain  is  proportional  to  the  applied  stress,  and  further  assume  that  materials  are  isotropic  and 
nearly  incompressible  in  nature.  The  description  of  the  constitutive  relationship  between  stress 
and  strain  is  ultimately  expressed  in  terms  of  the  elasticity  parameters  E  (Young’s  modulus) 
and  v  (Poisson’s  ratio). 

A  finite  element  representation  of  the  model  is  constructed  from  the  source  image. 
Elements  of  the  mesh  are  grouped  using  a  K-means  algorithm  by  initializing  a  number  (N)  of 
seed  points  that  are  the  centers  of  the  clusters  and  iteratively  minimizing  their  summed  distance 
to  all  element  centroids  in  the  mesh.  This  process  defines  a  set  of  nearly  equally  sized  but 
spatially  non-uniform  regions  that  are  homogeneous  with  respect  to  their  material  properties 
and  establish  the  ‘resolution’  of  the  reconstructed  elasticity  image.  After  assigning  appropriate 
boundary  conditions  based  on  estimated  displacement  or  stress,  the  standard  Galerkin  method 
of  weighted  residuals  (Lapidus  and  Pinder  1982)  is  used  to  construct  a  matrix  system.  The 
solution  of  that  system  yields  displacements  that  are  used  to  deform  the  source  image.  A 
second  discretization  is  performed  by  binning  the  target  image  into  M  groups  of  contiguous 
voxels  termed  zones.  The  model-deformed  image  is  then  compared  to  the  target  by  summing 
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the  similarity  metric  evaluated  for  all  zones.  The  correlation  coefficient  (Fitzpatrick  et  al 
2000)  is  used  throughout  this  work  as  it  has  empirically  demonstrated  better  performance 
for  our  method  over  other  intensity-based  metrics  such  as  the  sum  of  squared  differences 
and  normalized  mutual  information.  Optimization  of  the  elasticity  parameters  is  taken  as  the 
minimization  of  the  objective  function: 

^  =  I^true  —  SestI2  (2) 


where  Strue  is  the  set  of  similarity  values  achieved  when  comparing  the  target  image  to  itself, 
iSest  is  the  similarity  between  the  target  and  model-deformed  source  images  using  current 
estimates  of  the  elastic  modulus  distribution  and  |  •  |  denotes  the  vector  L2-norm.  Note  that  by 
definition,  Strue  for  the  correlation  coefficient  has  a  constant  value  of  1 .  Differentiating  (2) 
with  respect  to  the  elasticity  distribution  and  setting  the  resulting  expression  equal  to  zero 
generates  a  series  of  nonlinear  equations  that  can  be  solved  using  the  Levenberg-Marquardt 
method: 


[JrJ  +  a/]{A£}  =  [J7']{5True  -  SEST] 


a  =  X 


N 

i— 1 


(3) 


where  J  is  the  Jacobian  matrix  of  size  M  x  N  and  A E  is  the  vector  of  updates  to  the  material 
property  distribution  defined  by  the  regions.  The  regularization  parameter  a  uses  an  empirical 
scalar  factor  X  as  determined  by  the  methods  described  in  Joachimowicz  et  al  (1991).  Each 
column  of  the  Jacobian  matrix  is  a  finite  difference  approximation  of  the  change  in  image 
similarity  over  all  zones  due  to  the  perturbation  of  a  single  material  property  region,  such  that 
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Modulus  values  contained  in  E  are  updated  by  A E  until  an  error  tolerance  on  the  relative 
objective  function  error  evaluation  is  reached  or  a  maximum  number  of  iterations  are 
completed.  Spatial  averaging  of  elasticity  values  in  the  model  and  solution  relaxation  between 
iterations  are  also  utilized  to  improve  the  stability  of  the  optimization. 


2.2.  Parallel  computing  framework 

The  transition  of  this  method  from  2D  to  3D  entails  a  much  higher  computational  overhead 
that  affects  all  parts  of  the  reconstruction.  The  mesh  needed  to  describe  the  entire  breast  as 
opposed  to  a  single  slice  is  at  least  20-40  times  greater  in  the  number  of  structural  components 
(nodes  and  elements),  and  the  model  must  account  for  an  additional  degree  of  freedom.  The 
resulting  system  of  equations  to  be  solved  is  thus  nearly  two  orders  of  magnitude  larger. 
The  finite  difference  approximation  of  each  column  of  the  Jacobian  matrix  requires  a  ‘forward 
solve’  consisting  the  biomechanical  model,  image  deformation  and  evaluation  of  the  similarity 
metric.  Because  this  must  be  done  for  every  elasticity  region,  attempting  to  adequately  sample 
the  spatial  domain  makes  the  building  of  this  matrix  the  primary  expenditure  of  computing 
resources. 
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In  order  to  achieve  a  reasonable  level  of  performance,  the  Message  Passing  Interface 
standard  for  parallel  processing  is  used  to  distribute  formation  of  the  Jacobian  among  a 
number  of  communicating  nodes  controlled  within  a  static  single  process,  multiple  data 
(SPMD)  scheme.  The  Portable  Extensible  Toolkit  for  Scientific  Computation  (PETSc)  (Balay 
et  al  2004,  1997)  has  provided  the  necessary  coding  base  for  interfacing  sparse  matrix 
system  solvers  with  our  C/C++  Gauss-Newton  optimization  routine.  This  design  scales 
readily  to  the  number  of  processors  available;  it  has  been  tested  on  a  homogeneous  cluster  of 
18  processors  (2.0  GHz  Pentium4  Xeon,  1  GB  RAM)  located  in  the  laboratory,  as  well  as  a 
heterogeneous  cluster  of  hundreds  of  processors  available  through  the  Vanderbilt  Advanced 
Computing  Center  for  Research  and  Education  project.  The  use  of  many  processors  is  capable 
of  producing  a  nearly  linear  speedup  and  otherwise  agrees  in  principle  with  the  performance 
impact  suggested  by  Ahmdahl’s  law  (Ahmdahl  1967). 


2.3.  Simulation  experiment  setup 

For  this  work,  a  simulation  experiment  is  defined  by  the  creation  of  an  idealized  target 
image  volume  from  a  deformation  achieved  by  the  specification  of  boundary  conditions  at  the 
surface  of  the  breast.  This  ensures  data  fidelity  in  order  to  effectively  evaluate  reconstruction 
performance  in  the  optimization  and  model  characteristics.  Two  image  volumes  of  human 
breast  were  made  available  to  further  test  the  modality  independence  of  the  algorithm.  The 
first  was  obtained  from  a  dedicated  breast  CT  scanner  (256  x  256  x  130,  voxel  size  0.6  mm3) 
as  described  in  Boone  et  al  (2006),  (2001)  and  Boone  and  Lindfors  (2006)  and  the  second 
from  a  Philips  Achieva  3.0  T  MR  unit  (256  x  256  x  98,  voxel  size  1.0  mm3)  using  a  clinically 
approved  transmit-receive  double-breast  coil  to  acquire  a  3D  T1 -weighted  exam  with  a  fat¬ 
nulling  inversion  pulse  (TR/TE/a/NEX  =  4.6  ms/2.3  ms/10°/l)  (Yankeelov  et  al  2007). 
The  surfaces  of  the  breast  were  segmented  (ANALYZE  6.0,  Mayo  Clinic,  Rochester,  MN)  to 
create  tetrahedral  meshes  composed  of  39  013  nodes  connected  in  214  163  elements  for  the 
CT  volume  and  20  623  nodes  and  111  142  elements  for  the  MR  volume.  A  2  cm  spherical 
tumor  was  synthetically  implanted  in  the  center  of  each  mesh  by  assigning  a  stiff  modulus 
to  appropriate  member  elements  that  was  six  times  higher  than  the  surrounding  material 
(Krouskop  et  al  1998,  Samani  et  al  2007).  Tissue  deformation  was  performed  by  creating 
a  set  of  displacements  calculated  to  approximate  a  Gaussian  stress  distribution  applied  to  a 
rectangular  area  on  the  lateral  surface  of  the  breast.  The  displacements  were  then  applied  to 
the  original  volumes  in  order  to  create  the  desired  target  images.  Figure  2  illustrates  the  setup 
of  the  simulation  data. 


2.4.  Reconstruction  experiments 

Reconstructions  using  spatial  a  priori  knowledge  of  the  location  and  size  of  the  inclusion 
were  first  performed  in  order  to  constrain  the  problem,  as  well  as  the  computational  expense 
of  the  Jacobian  matrix,  to  a  two-material  discrimination  of  relative  stiffness  (elastic  contrast). 
A  second  set  of  experiments  was  then  used  to  address  the  viability  of  the  method  to  perform 
a  generalized  detection  of  the  lesion  with  no  knowledge  of  the  actual  structure  of  the  domain. 
To  run  these  naive  reconstructions  for  the  CT  data  set,  3180  material  regions  and  733  voxel 
similarity  zones  were  partitioned,  while  in  the  MR  data  set,  3166  regions  and  768  zones 
were  used.  In  all  cases,  the  reconstruction  was  initialized  with  a  homogeneous  elasticity 
distribution,  and  the  value  of  Poisson’s  ratio  held  constant  at  v  =  0.485  to  represent  a  nearly 
incompressible  material. 
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Figure  2.  (a)  CT  data  set  and  (b)  MR  data  set  used  for  3D  MIE  simulations.  Surface  renderings 
of  the  image  volumes  (top  row)  and  meshes  (bottom  row)  are  shown  for  the  pre-  ( source )  and 
post-deformation  ( target )  scenarios. 


2.5.  Evaluating  boundary  condition  influence 

In  addition  to  image  acquisition,  the  other  major  input  to  the  reconstruction  algorithm  is  the 
delineation  of  boundary  conditions  on  the  region  of  interest  over  which  the  model  is  applied. 
While  relatively  easy  to  control  in  simulation,  in  a  real  clinical  situation,  this  presents  the 
challenge  of  accurately  determining  point  correspondences  between  the  source  and  target 
breast  surfaces.  The  effect  of  any  inaccuracies  is  cumulative,  as  errors  are  propagated  from  the 
model  to  the  image  deformation  and  finally  the  similarity  measurements.  In  previous  2D  work, 
manual  delineation  of  boundary  conditions  was  possible  with  guidance  and  correction  using 
standard  computer  input  devices  (i.e.  a  mouse).  However,  the  increased  complexity  of  mesh 
geometry  in  3D  necessitates  a  more  automated  technique  of  determining  correspondence 
between  two  surfaces.  Potentially  non- trivial  random  and/or  operator-dependent  noise  is 
introduced  into  any  generated  boundary  conditions.  Therefore,  the  following  experiments 
were  performed  to  examine  the  ability  of  the  algorithm  to  tolerate  various  types  of  mis- 
mappings. 


2.5.7.  Robustness  to  randomized  boundary  errors.  The  gold  standard  boundary  conditions 
used  to  create  the  simulated  target  image  volumes  were  deliberately  disrupted  to  examine  the 
effect  of  random  noise  on  reconstruction  fidelity.  A  series  of  magnitudes  ranging  from  0.01  to 
2.0  voxel  units  (mesh  coordinates  normalized  by  their  respective  spacing  in  image  space)  were 
applied  to  the  CT  and  MR  data  sets.  Therefore,  every  boundary  position  is  displaced  by  the 
same  amount  but  in  a  completely  unpredicted  manner,  as  illustrated  in  figure  3.  These  altered 
boundary  conditions  sets  were  utilized  in  the  reconstruction  of  the  a  priori  two-material  test 
case,  and  the  tolerance  of  the  method  was  evaluated  by  calculating  the  average  reconstructed 
elasticity  contrast  ratios  over  four  trials  of  each  level  of  noise,  with  deviations  less  than  20% 
from  the  true  stiffness  being  deemed  acceptable. 


2.5.2.  Feasibility  of  automated  boundary  condition  generating  methods.  Three  methods  of 
surface  registration  and  point  correspondence  were  considered  as  the  basis  of  a  semi-automated 
method  for  determining  boundary  conditions  input  to  the  reconstruction  algorithm.  Two  were 
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(a) 


(b) 


Figure  3.  Examples  of  distortion  due  to  additive  randomized  error.  For  effect,  noise  of  2.0  voxel 
units  is  shown  as  applied  to  the  gold  standard  boundary  condition  set  for  CT  (a)  and  MR  (b).  At 
these  extreme  levels,  the  smooth  surface  of  the  breast  as  originally  captured  in  figure  2  is  completely 
lost,  and  the  forced  reconfiguration  of  internal  elements  in  the  finite  element  mesh  adversely  affects 
all  aspects  of  the  reconstruction. 


specifically  developed  for  this  work  by  attempting  to  use  potential  energy  distributions  derived 
from  classic  PDEs  for  surface  matching,  and  the  other  is  a  free-form  geometrical  warping. 

If  the  flow  of  a  hypothetical  substance  over  both  the  source  and  target  breast  surfaces  is 
taken  to  be  a  conserved  process  and  modeled  using  potential  theory,  correspondence  can  be 
assigned  by  matching  areas  of  similar  energy  deposition,  that  is,  the  equivalent  level  sets.  The 
algorithm  for  the  PDE-based  surface  matching  methods  can  be  summarized  in  the  following 
steps: 

(1)  Determine  an  energy  distribution  for  each  surface.  Laplace’s  equation  is  commonly  used 
to  describe  the  steady-state  distribution  of  potential  energy  O  in  a  system: 


-V2d>  =  0. 


(5) 


Similarly,  the  diffusion  equation  describes  the  temporal  change  in  potential  over  a  region: 


(6) 


where  D  is  the  diffusion  coefficient.  Each  PDE  is  solved  over  both  breast  surfaces  (source 
and  target).  For  both  equation  (5)  and  (6),  the  nipple  is  assigned  as  an  area  of  high 
potential  energy.  Additionally,  with  equation  (5),  nodes  at  the  chest  wall  are  assigned  a 
value  of  0  in  order  to  obtain  a  non-trivial  solution,  whereas  the  propagating  front  produced 
by  (6)  is  artificially  halted  at  the  chest  wall  boundary.  While  both  PDE  solutions  similarly 
establish  an  energy  gradient  over  the  breast  surfaces,  their  application  in  the  following 
steps  results  in  more  apparent  differences. 

(2)  Determine  correspondence  between  energy  distributions.  From  the  solution  of  the  PDEs 
on  the  source  surface,  a  series  of  spatially  distributed  isocontours  representing  distinct 
potentials  are  determined.  For  each  level  set,  an  isocontour  of  equivalent  potential  energy 
is  found  on  the  target  surface  and  the  two  curves  matched  according  to  the  symmetric 
closest  point  method  described  by  Papademetris  et  al  (2002). 

(3)  Generate  boundary  conditions.  By  extracting  a  number  of  isocontours  of  different  values, 
the  resulting  point  correspondence  vectors  define  a  relatively  dense  3D  displacement  field. 
The  displacement  for  each  boundary  node  can  then  be  interpolated  from  the  set  of  its 
nearest  neighbors. 
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The  final  method  employs  thin-plate  splines  (Goshtasby  1988)  to  generate  a  set  of 
boundary  conditions.  In  this  well-established  method  of  non-rigid  transformation,  a  number 
of  control  points  with  known  correspondences  establish  constraints  on  the  deflection  of 
a  hypothetical  thin  sheet  of  material  in  order  to  best  warp  the  two  surfaces  together. 
Displacements  at  each  boundary  node  are  then  simply  interpolated  from  the  calculated  fit. 
For  these  simulation  experiments,  a  subset  of  boundary  nodes  was  used  to  represent  physical 
markers  on  the  breast  surface.  Forty  points  were  uniformly  distributed  over  the  CT  mesh 
and  eighty  for  the  MR  mesh  in  order  to  handle  the  more  highly  variegated  shape  of  the  latter 
data  set. 

The  automated  methods  were  initially  evaluated  according  to  their  target  registration  error 
(TRE),  which  was  calculated  as  the  average  Euclidean  distance  between  the  generated  and  true 
boundary  conditions.  Because  the  deployment  of  these  fits  represents  a  more  correlated  form 
of  noise,  these  boundary  conditions  were  also  applied  to  the  two-material  (< a  priori )  scenario 
and  the  reconstructed  elasticity  contrast  values  compared  to  the  trials  of  additive  randomized 
error  for  which  the  magnitude  was  approximately  equal  to  the  TRE.  Finally,  a  mapping  of  the 
objective  function  space  was  performed  by  calculating  the  similarity  values  for  model-based 
image  deformations  over  a  range  of  elasticity  contrasts  from  0.5:1  to  30:1.  An  interpolating 
curve  was  fit  to  extract  the  minimum  objective  function  value  and  associated  contrast  ratio 
to  determine  a  theoretical  optimal  reconstruction  as  constrained  by  the  estimated  boundary 
conditions. 


3.  Results 


3.1.  MIE  reconstructions 


Because  the  use  of  a  priori  spatial  information  about  the  inclusion  limits  the  reconstruction 
to  a  two-material  system,  the  fidelity  of  the  reconstruction  is  simply  evaluated  by  examining 
the  elastic  contrast  between  the  inclusion  and  the  normal  tissue  of  the  breast  (ideal  of  6:1). 
Figure  4  demonstrates  the  behavior  of  the  algorithm  in  optimizing  the  objective  function  while 
successfully  characterizing  the  stiffness  of  the  inclusion  to  within  5%  of  the  actual  value 
(6.02:1  and  6.21:1  for  the  CT  and  MR  data  sets,  respectively). 

The  fidelity  of  the  generalized  reconstruction  experiments  (using  no  a  priori  knowledge  of 
the  domain)  was  primarily  evaluated  on  its  ability  to  detect  the  presence  of  an  inclusion  based 
on  classification  of  the  material  property  distribution  as  well  as  the  retrospective  accuracy  of 
localizing  the  lesion.  The  final  elasticity  values  were  treated  as  a  Gaussian  mixture  of  two 
classes  and  separated  by  a  threshold  established  via  the  method  described  in  Otsu  (1979). 
The  likelihood  of  discriminating  a  lesion  in  the  resulting  elasticity  image  was  found  using  the 
contrast- to-noise  ratio  (CNR)  as  defined  by  Bilgen  (1999)  and  Doyley  et  al  (2003): 


CNR  = 


/20l  -  /xb)2 


(7) 


where  /./  and  a2  are  the  sample  mean  and  variance  of  a  material  property  distribution  and  the 
subscripts  L  and  B  denote  the  lesion  and  normal  material  types,  respectively.  As  a  quantitative 
assessment  of  the  localization  of  the  lesion,  the  positive  predictive  value  of  correctly  identifying 
a  lesion  material  within  the  known  segmented  region  of  the  inclusions  was  also  calculated  as  a 
‘quality  of  reconstruction  score’  (QRS).  This  value  is  significant  because  identification  of  the 
lesion  border  and  material  classification  are  done  independently,  so  user  knowledge  of  the  test 
scenario  does  not  influence  the  performance  of  the  measure.  The  ‘true  positive’  (TP)  elements 
of  the  mesh  are  counted  as  the  number  correctly  identified  as  tumor  and  lying  within  the 
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Figure  4.  Optimization  behavior  of  reconstructions  using  a  priori  knowledge  of  the  inclusion 
location.  For  the  CT  simulation,  the  objective  function  evaluation  (normalized  to  the  initial 
dissimilarity  value  of  a  homogeneous  elasticity  distribution)  and  elastic  contrast  over  several 
iterations  of  the  algorithm  are  shown  in  panels  (a)  and  (b),  respectively.  Similarly,  this  behavior  for 
the  MR  data  set  is  displayed  in  (c)  and  (d).  In  each  case,  the  minimum  value  is  achieved  quickly 
and  stably,  with  the  corresponding  contrast  ratio  matching  the  true  value  of  6:1  very  closely  (6.02: 1 
and  6.21:1  for  CT  and  MR,  respectively). 


Table  1.  Evaluation  of  reconstruction  fidelity  for  lesion  detection. 


CNR 

QRS  (%) 

Max  CR  (x:l) 

Optimal  CR  (x:l) 

CT  3.55 

99.4 

2.66 

3.01 

MR  3.93 

99.7 

2.02 

2.26 

Max  CR:  maximum  elasticity  contrast  between  lesion  and  normal  tissue  in  naive  reconstruction. 
Optimal  CR:  optimal  elasticity  contrast  after  accounting  for  overlap  in  elasticity  region  partitioning. 


known  segmentation  of  the  lesion,  while  the  ‘false  positive’  (FP)  elements  are  those  identified 
as  tumor  but  in  an  incorrect  location.  Thus,  the  calculation  of  QRS  is  simply  TP/ (TP  +  FP). 
Cutoffs  for  successful  detection  and  localization  were  set  at  CNR  ^  2.2  as  noted  by  Doyley 
et  al  (2003)  and  QRS  >  80%  as  empirically  determined  by  a  prior  study  in  2D  MIE  work 
(Ou  et  al  2006a,  2006b),  and  both  the  CT  and  MR  reconstructions  successfully  identified  the 
embedded  lesions  according  to  these  criteria  (see  table  1). 

The  peak  modulus  contrast  value  of  a  reconstruction  was  calculated  by  taking  the  ratio  of 
the  average  elasticity  for  manually  selected  homogeneous  regions  of  the  approximately  equal 
area  known  to  be  representative  of  the  two  materials.  As  reported  in  table  1  and  visualized 
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Table  2.  Effect  of  applied  random  boundary  condition  noise  on  objective  function  space  and 
reconstructed  elasticity  contrast  ratio.  The  respective  ranges  where  a  cutoff  in  reconstruction 
tolerance  was  observed  are  listed  for  each  simulation  set. 


CT 

MR 

Randomized 

Mean  optimal 

Randomized 

Mean  optimal 

vector  magnitude 

elasticity  contrast 

vector  magnitude 

elasticity  contrast 

(voxel  units) 

value  (x:l) 

(voxel  units) 

value  (x:l) 

0.1 

5.62  ±  0.421 

0.01 

6.33  ±  0.096 

0.2 

5.70  ±  0.588 

0.02 

6.75  ±  0.058 

0.3 

5.97  ±  0.846 

0.03 

6.93  ±  0.634 

0.5 

2.36  ±  0.393 

0.05 

7.60  ±  0.821 

1.0 

2.47  ±  0.266 

0.1 

9.35  ±  1.27 

2.0 

2.17  ±  0.422 

0.2 

11.3  ±0.866 

Table  3.  Reconstruction  performance  as  affected  by  semi- automated  boundary  condition 
generation  methods.  The  mean  error  of  surface  registration  is  related  to  the  accuracy  of 
characterizing  the  lesion  stiffness. 


Method 

CT 

MR 

TRE 

(mm) 

Elasticity 
contrast  (x:l) 

TRE 

(mm) 

Elasticity 
contrast  (x:l) 

Diffusion 

1.5 

17.5 

0.61 

348 

Laplace 

0.52 

5.02 

0.48 

673 

Thin-plate  spline 

0.26 

5.66 

0.023 

6.26 

in  figures  5  and  6,  the  characterization  of  the  relative  stiffness  is  less  than  the  true  elasticity 
contrast  by  nearly  a  factor  of  3  in  both  cases.  This  reveals  a  difficulty  with  large  inverse 
problems  in  3D  where  a  need  for  reasonable  performance  can  lead  to  a  tradeoff  in  accuracy. 
By  choosing  approximately  3200  regions  to  cover  the  domain  of  the  breasts  for  the  naive 
reconstructions,  the  number  of  degrees  of  freedom  presented  to  the  optimization  scheme  is 
quite  high.  However,  this  is  also  relatively  coarse  in  the  sense  of  visualizing  the  reconstruction, 
as  it  roughly  corresponds  toal5  x  15  x  15  image  volume.  Because  the  elasticity  regions 
do  not  conform  perfectly  to  the  actual  lesion  borders  and  furthermore  comprise  both  the 
tumor  and  healthy  tissue,  it  seemed  reasonable  to  surmise  that  in  this  mis-estimation  of  spatial 
extents  the  algorithm  was  forced  to  attempt  a  best-fit  compromise.  To  test  this  hypothesis,  we 
agglomerated  all  regions  in  the  original  partitioning  that  overlapped  the  tumor  and  then  ran 
the  reconstruction  again  as  a  two-material  characterization.  Upon  inspection,  this  re-grouping 
was  clearly  a  larger  entity  than  the  tumor  itself  (closer  to  3  cm  in  diameter)  and  resulted  in  a 
shift  of  the  global  optimum  to  a  lower  elasticity  contrast.  In  effect,  the  model  reacted  to  this 
new,  oversized  tumor  by  reducing  its  stiffness  in  order  to  achieve  the  proper  image  similarity 
match.  When  viewed  in  light  of  this  analysis  as  summarized  in  table  1,  the  elasticity  contrast 
found  by  the  naive  reconstruction  is  then  actually  quite  accurate. 


3.2.  Evaluating  boundary  condition  influence 

3.2.1.  Robustness  to  randomized  noise.  Table  2  demonstrates  that  as  the  magnitude  of 
the  applied  randomized  noise  vectors  was  increased,  changes  in  the  reconstructed  elasticity 
contrast  reflected  a  decreased  ability  to  achieve  a  successful  result  (recall  that  the  correct  ratio 
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Figure  5.  Reconstruction  used  for  lesion  detection  in  the  CT  data  set.  (a)  Orthogonal  views 
taken  through  the  center  of  the  elasticity  image  volume  are  shown  along  with  a  projection  surface 
rendering  (lower  right).  The  simulated  inclusion  implanted  in  the  mesh  is  visually  distinguished 
from  the  surrounding  tissue.  The  color  bar  indicates  the  range  of  elasticity  values  (~7^t2  kPa) 
designated  by  the  reconstruction,  with  higher  (stiffer)  values  shown  in  the  white  end  of  the  grayscale 
mapping,  (b)  Transect  plots  through  the  center  of  the  volume  along  the  cardinal  directions  show 
the  profile  of  elasticity  contrast  (dotted  lines)  overlaid  by  the  true  profile  of  the  simulation  (solid 
lines). 


(a) 


(b) 


Figure  6.  Reconstruction  used  for  lesion  detection  in  the  MR  data  set.  (a)  Orthogonal  views 
taken  through  the  center  of  the  elasticity  image  volume  are  shown  along  with  a  projection  surface 
rendering  (lower  right).  Once  again,  the  inclusion  appears  to  have  a  recognizably  different  elasticity, 
with  values  on  the  color  bar  ranging  from  ~  10-57  kPa.  (b)  Transect  plots  through  the  center  of  the 
volume  along  the  cardinal  directions  show  the  profile  of  elasticity  contrast  (dotted  lines)  overlaid 
by  the  true  profile  of  the  simulation  (solid  lines). 


is  6:1).  For  the  CT  simulation,  on  average,  errors  of  0.5  voxel  units  or  greater  showed  a 
dramatically  reduced  ability  to  accurately  characterize  the  stiffness  of  the  lesion.  Similarly, 
though  at  a  much  smaller  scale,  the  MR  simulation  began  to  have  noticeable  difficulty  in 
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' - ■ - '  TRE  (mm) 
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Figure  7.  Three  candidate  automated  methods  for  MIE  boundary  condition  generation  applied  to 
simulation  CT  data.  Top  row,  from  left  to  right:  surface  deformations  calculated  from  diffusion 
energy  matching,  Laplace  solution  energy  and  thin-plate  spline  interpolation.  Bottom  row:  target 
registration  error  (TRE)  distribution  for  each  method  when  compared  against  the  gold  standard 
of  known  correspondence.  The  diffusion-based  mesh  is  both  qualitatively  and  quantitatively  the 
worst  performer.  The  Laplace  solution  appears  to  capture  the  shape  of  the  bladder  indentation 
more  precisely,  but  the  thin-plate  spline  has  the  best  overall  accuracy  in  determining  point 
correspondence. 

achieving  a  reasonable  reconstruction  at  noise  levels  of  0.05  voxel  units.  These  values  were 
taken  as  suitably  conservative  measures  for  evaluating  the  efficacy  of  boundary  conditions 
generated  by  the  semi- automated  methods. 

3.2.2.  Reconstruction  effects  of  generated  boundary  conditions.  The  accuracy  of  each 
automated  boundary  condition  technique  described  in  section  2.5.2  was  assessed  by  the  TRE 
with  respect  to  the  gold  standard  boundary  condition  set  and  its  ability  to  characterize  the 
elastic  contrast  in  the  two-material  reconstruction  test  case.  Figure  7  depicts  the  deformation 
fields  as  applied  to  the  CT  data.  Qualitatively,  the  displacements  found  by  the  diffusion 
method  are  quite  different  from  the  true  set,  while  the  results  from  the  solution  of  Laplace’s 
equation  and  the  thin-plate  spline  interpolation  appear  to  be  more  satisfactory.  As  shown  in 
table  3,  the  mean  TRE  of  the  three  methods  confirms  that  the  spline-based  method  has  the  best 
performance  (0.26  mm),  the  Laplace  method  next  (0.52  mm)  and  the  diffusion  method  being 
the  worst  (1.5  mm).  Inspection  of  figures  8  and  9  further  demonstrates  that  the  imposition 
of  an  inexact  boundary  condition  set  on  the  model  has  a  distinct  effect  on  the  optimization 
by  shifting  the  objective  function  minimum  value  to  a  different  optimal  elastic  contrast  ratio. 
Additionally,  the  convexity  of  the  objective  function  is  lost  in  the  cases  with  a  higher  TRE. 
The  differences  in  the  generated  boundary  condition  sets  for  the  MR  simulation  are  not  easily 
visualized  but  follow  a  similar  performance  trend  (TRE  of  spline  0.023  mm,  Laplace  method 
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Figure  8.  Mappings  of  objective  function  value  versus  elasticity  contrast  ratio  (tumor:breast) 
affected  by  the  boundary  condition  sets  generated  from  the  different  automated  methods  of  surface 
point  correspondence  as  applied  to  the  CT  data  set.  The  minimum  value  of  each  curve  corresponds 
to  the  altered  optimal  elasticity  contrast  when  constrained  by  the  inaccuracies  of  the  methods: 
(a)  diffusion,  (b)  Laplace  and  (c)  thin-plate  spline  interpolations.  The  ordinate  is  normalized  to  the 
initial  value  of  each  case.  The  global  minimum  of  (a)  is  out  of  range  of  the  plot. 

0.48  mm,  diffusion  0.61  mm).  For  both  simulations,  there  exists  a  direct  relationship  between 
a  low  TRE  and  increased  reconstruction  fidelity  in  characterization  of  the  elasticity  contrast 
of  the  lesion. 

4.  Discussion 

As  other  researchers  have  noted,  the  incorporation  of  a  priori  information  can  greatly 
enhance  the  performance  of  their  elastography  methods  (Doyley  et  al  2005,  2006).  We 
recognize  that  the  judicious  use  of  information  regarding  lesion  morphology  as  obtained 
from  conjunctive  imaging  studies  and  post-processing  would  potentially  aid  MIE  as  well, 
especially  in  reducing  the  number  of  search  parameters  and  improving  initialization  of  the 
algorithm.  The  reconstructions  using  a  priori  spatial  knowledge  of  the  inclusion  were 
initially  intended  to  simply  illustrate  that  the  objective  function  space  formed  by  using  an 
image  similarity  metric  was  smooth  and  readily  traversed  by  the  algorithm  in  a  manner 
expected  for  a  Gauss-Newton  optimization.  However,  they  also  provide  a  stark  contrast  to 
the  naive  lesion  detection  test  cases,  which  were  performed  to  evaluate  the  inverse  problem 
framework  and  demonstrate  its  ability  to  analyze  the  full  3D  domain  of  the  breast.  The  results 
of  the  generalized  reconstructions  are  very  encouraging  in  having  successfully  identified 
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Figure  9.  Mappings  of  objective  function  value  versus  elasticity  contrast  ratio  (tumor:breast) 
affected  by  the  boundary  condition  sets  generated  from  the  different  automated  methods  of  surface 
point  correspondence  as  applied  to  the  MR  data  set.  The  minimum  value  of  each  curve  corresponds 
to  the  altered  optimal  elasticity  contrast  when  constrained  by  the  inaccuracies  of  the  methods: 
(a)  diffusion,  (b)  Laplace  and  (c)  thin-plate  spline  interpolations.  Again,  the  ordinate  is  normalized 
to  the  initial  value  and  should  not  be  interpreted  as  an  equivalent  scale  for  each  case.  The  global 
minima  of  (a)  and  (b)  are  out  of  range  of  the  plot. 


and  localized  the  inclusions.  Although  the  discretizations  of  the  meshes  did  not  achieve 
particularly  accurate  material  characterizations,  the  optimal  elasticity  contrast  as  dictated  by 
the  available  objective  function  was  matched  in  each  case  to  within  12%.  The  observation  that 
mis-estimation  of  the  lesion  extent  altered  the  underlying  test  scenarios  suggests  that 
investigating  methods  of  dynamically  adjusting  region  assignment  could  facilitate  shape 
resolution  and  concomitantly  better  elasticity  contrast  ratio  values. 

In  translating  MIE  and  its  associated  technologies  to  a  clinical  setting,  a  number  of  factors 
must  be  considered  for  realistic  deployment.  From  an  implementation  and  performance 
perspective,  the  large  size  of  the  inverse  problem  necessitated  the  careful  selection  of  matrix 
solvers  and  programming  of  parallel  computing  routines  that  proved  effective  with  the 
availability  of  a  number  of  processors.  Initial  predictions  based  on  sequential  execution 
times  needed  to  handle  the  high  degrees  of  freedom  in  the  naive  reconstructions  were  thus 
reduced  from  two  weeks  to  several  hours.  Additional  challenges  were  eventually  overcome  in 
the  pre-processing  load  of  image  segmentation,  model  generation  and  partitioning  schemes. 

The  results  presented  in  this  paper  also  further  our  understanding  of  how  the  loss  of 
input  data  quality,  whether  through  design  limitations  or  unpredictable  factors,  could  have  a 
significant  impact  on  the  end  reconstruction.  In  particular,  the  proper  application  of  accurate 
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boundary  conditions  plays  a  critical  role  in  MIE  reconstruction  success.  This  is  due  to  the 
link  between  surface  shape  matching  and  subsequent  interpolation  of  internal  displacements 
in  affecting  sub-surface  image  intensities  and  similarity  measurements.  The  results  of  the 
boundary  condition  noise  experiment  are  interesting  because  they  indicate  that  some  level  of 
improper  localization  of  surface  point  correspondence  is  reasonably  tolerated  by  the  algorithm. 
However,  perturbations  greater  than  an  empirically  observed  threshold  can  impair  its  ability 
to  determine  the  underlying  elasticity  distribution.  This  is  a  similar  result  to  prior  work  done 
in  2D  systems  for  which  successful  reconstructions  correlated  to  boundary  condition  selection 
errors  limited  to  half  a  pixel  length  (Ou  et  al  2006b).  It  also  confirms  that  randomizing  the 
vectors  for  the  additive  noise  experiments  poses  a  considerable  challenge  to  the  algorithm 
because  of  the  introduction  of  grossly  non-physical  deformations  in  the  finite  element  mesh 
that  decrease  the  stability  of  the  numerical  model.  We  observed  that  the  threshold  for  the 
MR  simulation  was  an  order  of  magnitude  less  than  that  of  the  CT  set  and  initially  seemed 
to  require  an  unfeasible  level  of  accuracy,  as  well  as  quite  a  few  more  fiducials.  These  key 
differences  are  likely  related  to  image  resolution  (the  MR  volume  had  fewer  slices  and  a 
larger  voxel  spacing)  and  to  the  inherent  differences  in  soft-tissue  contrast  between  the  two 
modalities.  Both  issues  present  interesting  questions  that  will  be  explored  in  future  work. 

The  implausibility  of  performing  manual  selection  on  all  boundary  nodes  of  a  3D  mesh 
(there  were  6319  points  for  the  CT  and  5416  for  the  MR  set)  underscores  the  importance 
of  finding  an  automated  method  for  determining  point  correspondences.  In  general,  energy 
matching  from  the  solutions  of  the  diffusion  and  Laplace  equations  yielded  boundary  condition 
sets  that  were  inadequate  for  reconstructing  a  proper  elasticity  contrast.  This  can  be  partly 
explained  because  the  TRE  of  those  surface  registration  techniques  (as  compared  to  the  gold 
standard)  was  typically  greater  than  the  permissible  value  established  by  the  robustness  tests. 
The  primary  manifestation  of  these  poor  matches  was  that  the  model  often  had  difficulty  in 
obtaining  a  stable  solution.  Indeed,  only  the  boundary  conditions  generated  by  thin-plate  spline 
method,  which  had  the  least  error,  were  able  to  consistently  achieve  successful  reconstructions 
while  also  satisfying  the  putative  cutoffs.  Overall,  the  reconstruction  behavior  for  this  method 
was  consistent  to  within  6%  of  the  true  value.  This  appears  to  recommend  the  use  of  thin-plate 
spline  interpolation  as  a  strong  candidate  for  generating  boundary  conditions  for  MIE. 

5.  Conclusions 

In  this  work,  we  have  presented  the  first  fully  3D  realization  of  the  MIE  algorithm 
and  preliminary  evaluation  of  accompanying  strategies  for  automated  boundary  condition 
deployment.  The  use  of  parallel  processing  enabled  a  practical  implementation  of  a 
computational  problem  that  might  otherwise  prove  intractable.  Simulation  experiments 
demonstrate  the  viability  of  the  method  to  utilize  images  obtained  from  different  sources 
in  reconstructing  an  embedded  lesion  with  or  without  the  benefit  of  a  priori  information 
concerning  its  location  and  size.  We  have  also  characterized  the  robustness  of  the  elastography 
method  to  inaccuracies  in  boundary  condition  inputs  derived  from  either  random  noise 
or  by  surface  point  correspondence  methods.  These  results  should  prove  valuable  in  the 
customization  and  streamlining  of  data  acquisition  and  pre-processing  for  forthcoming  clinical 
tests. 
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ABSTRACT 

Multiple  skin  conditions  exist  which  involve  clinically  significant  changes  in  elastic  properties. 
Early  detection  of  such  changes  may  prove  critical  in  formulating  a  proper  treatment  plan.  However, 
most  diagnoses  still  rely  primarily  on  visual  inspection  followed  by  biopsy  for  histological  analysis.  As  a 
result,  there  would  be  considerable  clinical  benefit  if  a  noninvasive  technology  to  study  the  skin  were 
available.  The  primary  hypothesis  of  this  work  is  that  skin  elasticity  may  serve  as  an  important  method 
for  assisting  diagnosis  and  treatment.  Perhaps  the  most  apparent  application  would  be  for  the 
differentiation  of  skin  cancers,  which  are  a  growing  health  concern  in  the  United  States  as  total  annual 
cases  are  now  being  reported  in  the  millions  by  the  American  Cancer  Society.  In  this  paper,  we  use  our 
novel  modality  independent  elastography  (MIE)  method  to  perform  dermoscopic  skin  elasticity 
evaluation.  The  framework  involves  applying  a  lateral  stretching  to  the  skin  in  which  dermoscopic 
images  are  acquired  before  and  after  mechanical  excitation.  Once  collected,  an  iterative  elastographic 
reconstmction  method  is  used  to  generate  images  of  tissue  elastic  properties  and  is  based  on  a  two- 
dimensional  (2-D)  membrane  model  framework.  Simulation  studies  are  performed  that  show  the  effects 
of  three-dimensional  data,  varying  subdermal  tissue  thickness,  and  nonlinear  large  deformations  on  the 
framework.  In  addition,  a  preliminary  in  vivo  reconstruction  is  demonstrated.  The  results  are 
encouraging  and  indicate  good  localization  with  satisfactory  degrees  of  elastic  contrast  resolution. 

KEY  WORDS:  Elastography,  Dermoscopy,  Mechanical  Properties,  Finite  Element,  Image  Similarity, 
Elasticity  Imaging 


1.  INTRODUCTION 

Because  multiple  conditions  exist  which  involve  clinically  significant  changes  in  skin  elasticity, 
early  detection  of  such  changes  could  prove  critical  in  formulating  a  proper  treatment  plan.  However, 
most  diagnoses  still  rely  primarily  on  visual  inspection  followed  by  biopsy  of  suspect  areas  for 
histological  analysis.  Perhaps  the  most  apparent  application  would  be  for  the  differentiation  of  skin 
cancers,  which  are  a  growing  health  concern  in  the  United  States  as  total  annual  cases  are  now  being 
reported  in  the  millions  by  the  American  Cancer  Society1.  Of  the  three  major  types  of  skin  cancer 
reported  annually,  basal  cell  carcinoma  (BCC)  makes  up  approximately  800,000+,  squamous  cell 
carcinoma  (SCC)  cases  number  approximately  200,000+,  and  a  remaining  60,000+  melanoma  cases.  With 
respect  to  BCC,  approximately  5-10%  of  these  can  be  aggressive  and  infiltrate  the  surrounding  tissue  and 
sometimes  into  bone  and  cartilage.  It  rarely  metastasizes  but  can  cause  scars  and  disfigurement.  With 
respect  to  SCC,  early  detection  is  the  key  to  successful  treatment.  If  left  unchecked,  SCC  can  also  cause 
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disfigurement  and  typically  approximately  3-4%  of  cases  results  in  metastasis  which  is  usually  fatal. 
Melanoma  is  the  most  fatal.  Melanoma  is  malignant  and  if  left  unchecked,  it  will  spread  to  other  parts  of 
the  body,  becomes  difficult  to  treat,  and  can  be  fatal.  If  one  were  to  include  among  the  common  cancer 
statistics,  aggressive  BCC’s  and  metastatic  SCC’s,  skin  cancer  would  likely  be  the  second  most  prevalent 
among  newly  diagnosed  cancers. 

While  the  above  statistics  are  compelling,  when  speaking  to  lethality,  skin  cancer  is  less 
significant  than  its  more  fatal  counterparts  in  breast,  lung,  colon/rectal,  and  prostate.  However  when 
considering  the  economic  cost  that  skin  cancer  incurs  on  US  healthcare,  the  pursuit  of  skin  cancer 
characterization  has  considerable  merit.  The  ability  to  differentiate  benign  from  malignant,  and 
aggressive  from  non-aggresssive  skin  lesions  would  provide  considerable  savings  to  health  care  costs 
within  the  United  States.  Each  year  it  is  estimated  that  approximately  5-7  million  patients  undergo 
biopsies  for  skin  cancer  with  only  a  fraction  resulting  in  malignancies.  While  complete  multi-center 
biopsy  studies  have  not  been  performed,  one  study  that  took  place  documented  the  percentage  of  skin 
biopsy  specimens  that  were  malignant,  i.e.  they  termed  this  the  “malignancy  ratio”  2.  In  this  study,  the 
malignancy  was  approximately  40%  with  a  very  wide  variability  among  the  22  dermatologists  (13.4- 
86.6%)  that  participated  in  the  study.  If  one  considers  5  million  biopsies,  this  would  translate  to 
approximately  3  million  unnecessary  biopsies  per  year,  and  approximately  $300-600  million  in  saved 
expenditures  per  year  3’ 4.  Furthermore,  this  does  not  even  include  cosmetic  surgery  costs  in  the  case  of 
scarring  or  complications  associated  with  bleeding  and  infections.  If  an  inexpensive  imaging  device 
could  differentiate  lesions  with  reasonably  high  specificity  and  sensitivity,  it  would  have  considerable 
significance. 

It  should  also  be  noted  that  measuring  cutaneous  elasticity  is  also  potentially  important  to  medical 
areas  outside  of  clinical  dermatology.  For  instance,  in  a  recent  study  of  100  women  receiving  hormone 
replacement  therapy,  Pierard  et  al  demonstrated  a  positive  correlation  between  bone  mass  density  and 
skin  elasticity  5’  6.  Another  study  conducted  by  Yoon  et  al  showed  a  similar  relevance  for  patients 
afflicted  with  diabetes  mellitus  7.  Further  uses  for  evaluating  skin  elasticity  range  from  surgery 
(minimization  of  scarring,  chronic  graft  versus  host  disease)  to  rheumatology  (scleroderma,  lupus)  to 
obstetrics  (striae  development  in  pregnancy) 8'12. 

With  respect  to  diagnostic  technological  advances,  developments  have  been  concerned  with 
obtaining  a  better  view  of  the  skin,  either  via  improved  optics  (i.e.  dermoscopy)  or  by  more  advanced  and 
novel  imaging  systems  ranging  from  high-frequency  ultrasound  to  confocal  laser  microscopy  13, 14.  Other 
strategies  involving  electrical  impedance  mismatch  15,  Raman  spectroscopy  16,  and  cytological  smears  14 
have  also  been  forthcoming.  In  lieu  of  these  methods  that  capitalize  on  electrical,  optical,  and 
biochemical  phenomena,  we  have  chosen  to  pursue  an  approach  which  is  based  on  analyzing  mechanical 
behavior  of  the  skin.  Detecting  changes  in  tissue  by  palpation  and  then  associating  them  with  a  disease 
state  has  a  long-standing  history  in  clinical  medicine,  and  utilizing  changes  in  the  mechanical  properties 
to  specifically  characterize  the  skin  does  have  precedent.  A  thoughtful  review  by  Edwards  and  Marks 
discusses  the  complex  mechanical  behavior  of  skin  when  subjected  to  in  vitro  and  in  vivo  testing  17.  Their 
review  highlights  extensive  methodologies  being  used  to  quantify  skin  properties  (e.g.  uni-axial  and  bi¬ 
axial  extensometry,  torsion  stimulators,  indentometery,  ballistometric  tests,  shear  wave  application 
devices,  dynamic  suction  methods,  ultrasonics,  and  electrodynamometry)  and  also  emphasizes  the 
difficulties  in  comparing  across  these  methods.  One  of  the  authors'  conclusions  is  that  the  need  for 
quantitative  reproducible  methods  to  assess  skin  health  is  necessary  given  the  considerable  subjectivity  in 
clinical  analysis. 

Following  previous  work  in  18,  we  are  developing  a  new  method  we  have  termed  ‘  modality  - 
independent  elastography’  (MIE)  that  combines  image  processing  with  an  inverse  problem.  More 
specifically,  image  similarity  metrics  routinely  used  with  image  registration  methods  are  recast  to  satisfy 
an  inverse  elasticity  problem  framework  whereby  mechanical  properties  within  a  biomechanical  model  of 
deforming  tissue  become  the  driving  parameters  for  improved  image  similarity.  In  this  way,  MIE 
circumvents  two  potential  limitations  of  current  elastographic  techniques.  First,  it  is  not  inherently 
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dependent  on  pre-processing  steps  such  as  homologous  feature  selection  and  tracking  which  drive  active 
contour  models  '  and  other  traditional  displacement-based  iterative  methods  '  ,  although  such 
techniques  can  aid  in  the  determination  of  boundary  conditions.  Secondly,  because  it  is  an  image 
processing  methodology,  MIE  is  not  reliant  on  a  particular  imaging  modality  (such  as  in  ultrasound  and 
magnetic  resonance  elastography)  as  long  as  the  acquired  images  provide  sufficient  pattern  to  allow  for 
comparison.  Building  upon  recently  with  a  multi-resolution  implementation25,  this  paper  presents 
analysis  using  a  tissue  model  that  incorporates  geometric  nonlinearities,  the  effects  of  the  three- 
dimensional  nature  of  the  problem  which  include  varying  subcutaneous  layer  thicknesses,  and  varying 
lesion  depth.  In  addition,  a  relatively  crude  preliminary  in  vivo  result  is  also  demonstrated. 

2.  MATERIALS  AND  METHODS 
2.1.  Modality  Independent  Elastography 


Fig  1.  MIE  framework. 

MIE  begins  with  the  acquisition  of  an  image  series  consisting  of  an  image  acquired  prior  to  and 
after  an  applied  deformation.  The  method  is  “independent”  because  it  does  not  require  any  specific 
imaging  modality  but  rather  that  sufficient  pattern  is  present  within  the  acquired  images  such  that  material 
properties  can  be  assessed  from  pattern  changes  within  the  acquisition-pair.  This  is  a  similar  requirement 
for  many  methods  of  non-rigid  image  registration.  At  its  core  MIE  is  an  image  analysis  method  whereby 
a  model-generated  deformation  field  is  applied  to  the  pre-deformed  image  series  {source)  and  compared 
to  its  acquired  deformed  counterpart  {target).  This  comparison  is  nested  within  a  Levenberg-Marquardt 
optimization  routine  such  that  the  material  properties  become  the  parameters  of  interest  in  matching  the 
model-deformed  source  image  to  the  acquired  target  image.  The  methods  and  development  of  this 
technique  have  been  reported  in  detail  elsewhere18,25'29. 

With  respect  to  the  optimization  framework  for  MIE,  it  can  be  represented  as  a  least  squared  error 
objective  function: 

4>(e)=  min{|s(E,)- s(ee|2}  (1) 

where  s(Et )  is  the  similarity  value  achieved  when  comparing  the  target  image  to  itself  (i.e.  the  maximum 

value  for  the  similarity  measure,  so  for  the  correlation  coefficient  S^Et)=  1)  and  s(ee)  is  the  similarity 
between  the  model-deformed  source  image  and  the  acquired  target  image  using  the  current  estimate  of  the 
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elastic  modulus.  Eq.  (1)  is  optimized  by  setting  the  partial  derivative  equal  to  zero  and  solved  using  a 
Levenberg-Marquardt  approach: 


[[JT][J]  +  a[l]]^E}=  [jt^(§,)-S(ee)}  (2) 

dSiE 

is  the  M  x  N  Jacobian  matrix  of  the  form  J  = - W 

8E 

measurements  (i.e.  zones)  being  made  and  N  is  the  number  of  material  property  regions.  Because 
NT][J]  (an  approximation  to  the  Hessian  matrix)  tends  to  be  ill-conditioned,  it  is  regularized  with  an 
empirically  determined  a  parameter  found  in  the  standard  Levenberg-Marquardt  approach30.  The 
determination  of  this  regularization  parameter  is  described  in  31.  The  multi-resolution  framework  to  MIE 
is  shown  in  Fig.  1.  The  methodology  used  involves  a  hierarchical  material  parameter  resolution  series. 
This  has  been  reported  elsewhere  25, 27  and  has  been  shown  to  assist  in  avoiding  local  minima  that  are 
associated  with  the  decorrelation  of  patterned  data. 

2.2.  Model  for  Skin 


where  [J] 


^  and  M  is  the  number  of  similarity 


One  critical  component  within  all  model-based  inverse  problem  frameworks  is  the  selection  of 
the  computer  model  to  represent  the  continuum  of  interest.  In  previous  phantom,  simulation,  and  in  vivo 
studies,  we  have  elected  to  employ  a  plane  stress  linear  elastic  model  to  simulate  the  skin.  This  model  is  a 
two  dimensional  approximation  of  the  three  dimensional  system  which  assumes  a  symmetric,  isotropic, 
thin  specimen  in  equilibrium,  and  stresses  that  are  constrained  to  lie  within  the  plane.  These  assumptions 
simplify  Cauchy’s  law  from  36  stiffness  constants  to  2  and  use  the  equation: 

V  •  a  =  0  (3) 

where  a  is  the  2D  Cartesian  stress  tensor  and  is  defined  below  followed  by  the  plane  stress  constitutive 
equations: 
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where  E  is  Young’s  modulus,  and  v  is  Poisson’s  ratio.  For  this  work,  we  have  determined  that  a  Poisson’s 
ratio  of  0.485  for  our  skin  phantoms  and  tissue  work  has  performed  reasonably  well.  This  value  would 
correlate  with  a  -32:1  ratio  of  the  Lame'  constants  (k  :p,  with  pthe  shear  modulus,  and  X  being  the 
second  Lame'  constant)  which  is  reasonably  below  the  convention  for  Poisson  locking  (sometimes  called 
mesh  locking  and  typically  has  k  »  \i )  although  one  could  argue  that  hyperelastic  models  may  be  the 
more  appropriate  model  and  will  need  to  be  explored  in  the  future.  Using  the  Galerkin  weighted  residual 
method  to  integrate  and  solve  this  set  of  partial  differential  equations,  a  finite  element  framework  is 
generated  and  can  be  solved  to  represent  a  displacement  field  for  a  given  distribution  of  Young’s 
modulus. 

In  this  paper,  the  plane  stress  model  has  been  enhanced  to  incorporate  the  full  Green-Lagrange 
strain  tensor  as  defined  by: 


1  f  SUj 


Use  of  this  tensor  description  abandons 
compatible  with  large  deformations.  The 
theory  can  be  seen  in  the  2D  simulations 


SUj  3uk  3uk ) 

— -  +  — ^  .  (6) 
dx,  dx,  3xjy 

the  traditional  small-strain  approximation  in  favor  of  one 
difference  in  solutions  between  small  and  large  deformation 
shown  in  Fig.  2.  Fig.  2  compares  the  boundary  shape  of  a 
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compressive  strain. 


Hookean  linear  elastic  membrane  to  a  2D  Hookean  geometrically  nonlinear  elastic  model.  For  these 
simulations,  a  plane  stress  approximation  was  performed  for  comparison  only  of  the  solution  due  to 
nonlinear  terms  (i.e.  a  -30-40%  compressive  stress  could  never  be  applied  to  a  thin  specimen,  the 
material  would  buckle  well  before).  The  important  aspect  to  notice  is  that  more  necking  and  less  bulging 
occur  in  the  nonlinear  than  the  linear  model  in  tension,  and  compression,  respectively.  With  respect  to  the 
difference  in  the  linear  model  among  Fig.  2a  and  2b,  2b  is  the  reverse  of  2a  (this  is  a  characteristic  of 
linear  theory).  However,  the  lack  of  this  symmetry  for  the  nonlinear  model  is  characteristic  of  the  Green- 
Lagrange  strain  tensor  and  is  caused  by  the  interplay  between  the  linear  and  quadratic  terms  in  the  tensor. 
In  this  paper,  the  results  from  the  linear  and  nonlinear  model  will  be  compared. 


Fig  3.  Phantom  membrane  stretches  (a)  baseline,  (b)  1  cm  stretch,  (c)  2  cm  stretch. 


2.3.  Nonlinear  Model  Experiments 

In  previous  work,  a  phantom  of 
simulated  skin  was  generated  from  using  a 
polyurethane  membrane  25  cm  long,  15  cm 
wide,  and  approximately  2  mm  thick  containing 
two  types  of  materials.  The  bulk  material,  representing  “normal”  skin,  was  chosen  to  be  Evergreen  10 
(Smooth-On,  Easton,  PA),  while  the  stiffer  Evergreen  50  was  used  to  create  a  5  cm  circular  inclusion  of 
full  depth  and  was  placed  at  the  center  of  the  phantom.  The  visible  surface  of  the  phantom  was  modified 
by  drawing  a  regular  grid  of  horizontal  and  vertical  lines  spaced  about  1  cm  apart  in  either  direction.  Fig. 
3  shows  the  skin  phantom  as  acquired  by  an  optical  camera  at  the  baseline,  1  cm,  and  2cm  stretch  levels. 
Separate  mechanical  tests  on  the  membrane  were  conducted  using  the  EnduraTEC  ELF  3200  material 
tester  (EnduraTEC  Systems  Group,  Bose  Corporation,  Minnetonka,  MN) 25 .  Table  1  reports  the  expected 
Young’s  Modulus  contrast  ratio  based  on  a  Hookean  solid  fitted  to  each  respective  stretch/strain  level. 
The  scale  of  lesion  to  field-of-view  size  is  the  anticipated  aspect  ratio  for  a  dermoscopic  system.  This 


Stretch 

CR  (Linear) 

CR  (Nonlinear) 

1  (0.5  cm) 

5.7 

5.0 

2  (1.0  cm) 

5.0 

4.8 

3  (1.5  cm) 

4.6 

4.7 

4  (2.0  cm) 

4.2 

4.4 

Table  1.  Expected  elastic  contrast  ratio. 
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image  data  was  used  as  an  input  to  the  MIE 
framework,  and  reconstructions  comparing  the 
linear  and  nonlinear  model  were  generated. 

2.4.  Three-Dimensional  Model  Experiments 

In  an  effort  to  test  the  MIE  algorithm  with 
more  realistic  images  of  the  skin  as  would  be 
expected  in  the  clinic,  a  simulation  study  was 
performed  on  an  image  obtained  from  the 
Dermatology  Image  Atlas  project 
(www.dermatlas.org,  “melanoma_l_0405 1 0”, 
contributed  by  Eric  Ehrsam,  M.D.,  Fig.  4a). 

Dermoscopic  images  provide  the  challenge  of 
relying  on  the  natural  pattern  instead  of  the  structured  grid  used  in  the  membrane  phantom  experiments 
(although  the  skin  could  be  printed  upon  with  an  ink  grid).  Previous  simulation  work  on  this  image  can 
be  found  in 25 . 

One  possible  critique  of  this  dermoscopic  framework  is  the  fear  that  underlying  layers  would 
confound  the  approach.  As  a  result,  we  have  modified  the  simulation  study  reported  in  25  (e.g.  Fig.  4b)  to 
be  more  realistic.  Fig.  5a-5c  shows  the  setup  of  the  new  simulation  study.  In  this  study,  six  10  cm  x  10 
cm  domains  were  constructed  which  had  different  layered  structures.  Three  domains  were  6mm  in  total 
thickness  (-4mm  epidermis/dermis,  2 mm  subcutaneous)  and  had  1cm  central  inclusions  varying  in  depth 
1mm,  2mm,  and  3mm,  respectively.  The  remaining  three  domains  were  the  same  except  that  the 
subcutaneous  layer  thickness  was  increased  to  7mm.  The  selection  of  subcutaneous  layer  thickness  32, 33 
and  stiffness  values34  were  based  on  the  literature.  Using  this  domain,  a  mechanical  aperture  that 
stretched  the  skin  was  simulated,  3D  displacements  were  calculated  (Fig.  5b),  and  then  projected  on  to  the 


b 


mm 


D=1cm 


Dermis  1 - — 1.2,3  mm 

1  — 
▼  T 

Subcutaneous  layer 

Fig.  5.  (a)  Mesh  with  refinement  in  lesion  area,  (b)  simulated  compression  of  skin  by  device,  and  (c)  depth 
variation  (transect  T  is  shown  in  (a)). 


original  2D  mesh.  A  new  set  of  simulated  images  (deforming  Fig.  4)  was  generated  and  then  used  as 
input  to  the  2D  MIE  algorithm. 

2.5  Noninvasive  In  Vivo  Experiments 


In  addition  to  the  simulation  experiments,  a  second  camera  and  deformation  system  was 
constructed  using  a  Sony  XCD-X710CR  CCD  camera  with  a  Schneider  12  mm  lens  and  a  circular 
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polarizing  filter  for  the  optics  mounted  over  a  crude  spring-loaded  standard  set  of  pincers  that  were 
pressed  against  the  skin  surface  and  bound  with  commercial  adhesive.  Digital  images  (800  x  600 
resolution)  of  a  common  nevus  of  palpable  stiffness  that  was  2-3  mm  in  diameter  and  located  on  the  volar 
aspect  of  the  left  forearm  of  a  male  volunteer  were  acquired  in  both  relaxed  and  stretched  states.  It  was 
determined  that  the  reflection  and  scattering  of  ambient  fluorescent  lighting  interfered  with  this  particular 
setup  and  affected  the  input  image  quality,  so  an  artificial  grid  was  imposed  to  counteract  specularity. 


Stretch  1 


Stretch  2 


4 


3 


Fig  6.  (a)  Linear  Hookean  reconstmctions  for  each  stretch  state,  and  (b)  the  geometrically  nonlinear  Hookean 
reconstmctions.  Contour  shows  inclusion  border. 


Images  were  acquired  in  the  baseline  and 
stretched  state  and  given  to  the  MIE  algorithm. 

3.  RESULTS 

Fig.  6  illustrates  the  reconstructions  of 
the  phantom  membrane  system  using  the  linear 
and  nonlinear  models,  respectively.  Fig.  7  reports 
how  the  objective  function  varies  between  the 
linear  and  nonlinear  model  reconstructions  of  Fig. 
6.  Fig.  8  demonstrates  the  dependence  of 
resolving  stiffness  (3:1,  6:1,  12:1)  on  depth  (1,2, 

Fig  7.  Objective  function  difference  between  linear  and  ^  mm)  f°r  the  18  combinations  (combinations  of 
nonlinear  models  over  all  stretches  and  resolutions.  ^  depths,  3  contrast  ratios,  and  2  subcutaneous 

layer  thicknesses).  Table  2  reports  an 
approximate  contrast  ratio  of  the  lesion-to-bulk  material  for  each  of  the  cases  shown  in  Fig.  8.  Fig  9. 
illustrates  the  results  from  the  in  vivo  experiment  conducted.  Fig.  9  shows  (a)  the  nevus  ,  (b)  the  finite 
element  grid,  and  (c,d)  reconstructed  elasticity  images.  Fig.  9c-d  illustrate  the  effect  of  incorporating 
increasing  degrees  of  a  priori  knowledge  of  the  actual  location  and  elasticity  distribution  into  the 
algorithm.  Fig.  9c  is  a  general  elasticity  reconstruction  of  the  nevus  with  lesion-related  initial  guess,  and 
Fig.  9d  maintains  the  spatial  prior  of  Fig.  9c  for  the  entire  duration  of  the  reconstruction  process. 


Strain  (%)  SO-  #  of  regions 


4.  DISCUSSION 


With  respect  to  Fig.  6,  the  most  important  features  to  note  are  a  very  subtle  improvement  in 
satisfying  the  inclusion  contour,  and  the  decrease  in  variability  of  the  contrast  in  the  geometric  nonlinear 
reconstructions.  This  is  consistent  with  the  behavior  in  Table  1  and  indicates  that  the  geometric  nonlinear 
reconstruction  is  performing  as  predicted.  One  thing  to  note  is  that  if  the  membrane  were  a  true  Hookean 
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2  mm  subcutaneous 


7  mm  subcutaneous 
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Fig.  8.  3-D  effects  on  elasticity  images  with  varying  lesion  depth  and  subcutaneous  layer  depth.  Ratio  on  left 
shows  actual  contrast  while  colorbar  shows  reconstructed.  Fig.  4  shows  the  lesion  contour. 


2mm 

SQ 

1  mm 

2  mm 

3  mm 

3:1 

2.12 

2.55 

2.79 

6:1 

3.08 

3.76 

4.00 

12:1 

4.11 

4.71 

5.15 

7  mm 

SQ 

1  mm 

2  mm 

3  mm 

3:1 

1.93 

2.27 

2.34 

6:1 

2.87 

3.49 

3.55 

12:1 

4.03 

4.28 

4.53 

Table  2.  Lesion-to-bulk  contrast  ratio  at 
each  lesion  depth  (1mm,  2mm,  3  mm) 
and  target  contrast  level  (3:1,  6:1,  12:1) 
for  2mm  (top)  and  7mm  (bottom) 
subcutaneous  thicknesses. 

nonlinear  solid,  the  CR  would  be  the 
same  across  all  stretches,  i.e.  strain 
levels.  While  the  variability  across 
levels  was  reduced  between  the 
linear  and  nonlinear  models  for 
strain,  the  data  supports  a  more 
complex  constitutive  model  for  this 
phantom  system.  Fig.  7  also  demonstrates  some  consistency  with  the  extension  to  a  nonlinear  model. 
Here  we  see  that  at  the  largest  stretch  values,  the  difference  between  the  two  models  is  highest  with  the 
nonlinear  model  outperforming  the  linear. 

With  respect  to  Fig.  8,  and  Table  2,  despite  not  reaching  the  expected  contrast  ratio  values,  each  result 
localized  the  lesion  and  demonstrated  sensitivity  to  depth  and  contrast.  With  respect  to  contrast  ratio,  this 


Fig.  9.  (a)  Dermoscopic  scale  interrogation  (magenta  border  is  about 
3  mm  in  diameter),  (b)  reconstmction  mesh,  (c)  general  elasticity 
reconstruction,  (d)  incorporation  of  lesion  border  as  a  priori 
information. 
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simulation  was  more  challenging  than  previous  work.  The  results  indicate  that  varying  subcutaneous 
thicknesses  has  a  very  modest  effect  in  confounding  lesion  reconstructions  for  the  contrast  regime  shown. 
Fig.  8  does  suggest  that  at  some  contrast  levels  there  may  be  overlap  (e.g.  a  3mm,  3:1  lesion  with  2mm 
subcutaneous  layer  may  be  confused  with  a  1mm,  6:1  lesion  with  7mm  subcutaneous  layer).  This  could 
be  potentially  confounding  if  left  in  the  absence  of  other  information.  However,  Horejsi  et  al.  and  Moller 
et  al.  32, 33  do  provide  guidance  for  estimating  the  subcutaneous  tissue  thickness  in  a  general  population. 
Additionally,  they  use  a  relatively  inexpensive  optical  device  to  make  these  measurements.  There  is  little 
doubt  that  using  this  a  priori  information  would  reduce  variability  in  interpretation.  It  is  interesting  to 
speculate  that  if  a  particular  lesion  type  was  of  a  known  stiffness  a  priori  then  perhaps  the  reconstruction 
contrast  could  be  used  to  estimate  lesion  depth.  Looking  across  the  grids  of  a  particular  property  contrast, 
there  is  a  definitive  change  in  the  reconstruction  as  a  function  of  depth. 

Fig.  9  demonstrates  a  preliminary  attempt  of  applying  the  framework  to  an  in  vivo  system. 
Considering  there  was  very  little  control  over  this  system  and  no  ground  truth  was  available,  these  results 
qualitatively  confirm  the  potential  utility  of  MIE  in  evaluating  an  area  of  skin  with  a  region  of  differing 
elasticity.  While  the  results  in  Fig.  9  c,d  are  not  completely  satisfying  at  this  time,  it  should  be 
emphasized  that  the  camera-system  employed  a  coarse  resolution,  and  the  skin-to-device  coupling  system 
was  extremely  crude.  Given  this  somewhat  ad-hoc  experiment,  the  detection  of  any  anomaly  in  the 
correct  region  is  encouraging  for  this  approach.  It  was  also  interesting  to  note  the  increase  in  performance 
by  adding  the  lesion  extents  to  the  information  provided  to  the  algorithm. 

5.  CONCLUSIONS 

The  results  from  the  experiments  above  demonstrate  an  ability  to  use  the  MIE  framework  within 
the  dermoscopy  setting.  Methodological  concerns  regarding  the  use  of  geometric  nonlinear  models, 
three-dimensional  effects,  and  in  vivo  conditions  are  addressed.  The  results  indicate  that  while  geometric 
nonlinearities  do  modestly  affect  reconstructions,  nonlinear  material  models  may  be  needed  to  correct  for 
remaining  discrepancies.  A  modest  improvement  was  shown  using  the  geometric  nonlinear  model 
especially  at  large  strain  values  which  is  consistent  with  the  theoretical  development.  The  three- 
dimensional  effects  of  lesion  depth  and  varying  subcutaneous  layer  thicknesses  are  assessed.  Lesion 
depth  does  affect  contrast  ratio  whereas  subcutaneous  layers  affect  the  reconstructions  to  a  significantly 
lesser  degree.  Avenues  to  detect  lesion  depth  in  the  presence  of  a  prior  knowledge  of  a  lesion’s  stiffness 
may  have  surgical  implications.  Preliminary  in  vivo  work  suggests  that  lesion  characterization  is  possible 
although  specificity  and  sensitivity  of  the  method  await  further  study  and  will  need  a  considerably  more 
robust  acquisition  system. 


6.  ACKNOWLEDGEMENTS 

This  work  was  supported  by  the  Whitaker  Foundation’s  Young  Investigator  Program  and  the  Whitaker 
Foundation’s  Continuation  Program.  This  work  was  also  supported  by  the  Congressionally  Directed 
Medical  Research  Program-Breast  Cancer  Research  Program. 

7.  REFERENCES 

1.  American  Cancer  Society,  ’’Cancer  Facts  and  Figures  2006,”  A.  C.  Society,  ed.  (Atlanta,  GA, 
2006). 

2.  A.  R.  Green,  G.  W.  Elgart,  F.  Ma,  D.  G.  Federman,  and  R.  S.  Kirsner,  ’’Documenting 
dermatology  practice:  Ratio  of  cutaneous  tumors  biopsied  that  are  malignant,”  Dermatol.  Surg.  30,  1208- 
1209  (2004). 


34 


3.  G.  J.  Chen,  C.  B.  Yelverton,  S.  S.  Polisetty,  T.  S.  Housman,  P.  M.  Williford,  H.  V.  Teuschler,  and 
S.  R.  Feldman,  'Treatment  patterns  and  cost  of  nonmelanoma  skin  cancer  management,"  Dermatol.  Surg. 
32,  1266-1271  (2006). 

4.  K.  A.  Freedberg,  A.  C.  Geller,  D.  R.  Miller,  R.  A.  Lew,  and  H.  K.  Koh,  "Screening  for  malignant 
melanoma:  A  cost-effectiveness  analysis,"  Journal  of  the  American  Academy  of  Dermatology  41,  738- 
745  (1999). 

5.  G.  E.  Pierard,  C.  Pierard-Franchimont,  S.  Vanderplaetsen,  N.  Franchimont,  U.  Gaspard,  and  M. 
Malaise,  "Relationship  between  bone  mass  density  and  tensile  strength  of  the  skin  in  women,"  Eur.  J. 
Clin.  Invest.  31,  731-735  (2001). 

6.  G.  E.  Pierard,  S.  Vanderplaetsen,  and  C.  Pierard-Franchimont,  "Comparative  effect  of  hormone 
replacement  therapy  on  bone  mass  density  and  skin  tensile  properties,"  Maturitas  40,  221-227  (2001). 

7.  H.  S.  Yoon,  S.  H.  Baik,  and  C.  H.  Oh,  "Quantitative  measurement  of  desquamation  and  skin 
elasticity  in  diabetic  patients,"  Skin  Res.  Technol.  8,  250-254  (2002). 

8.  C.  Braham,  D.  Betea,  C.  Pierard-Franchimont,  A.  Beckers,  and  G.  E.  Pierard,  "Skin  tensile 
properties  in  patients  treated  for  acromegaly,"  Dermatology  204,  325-329  (2002). 

9.  D.  Carrino,  J.  Sorrell,  and  A.  Caplan,  "Age-related  changes  in  the  proteoglycans  of  human  skin," 
Archives  of  Biochemistry  and  Biophysics  373,  91-101  (2000). 

10.  J.  A.  Clark,  J.  C.  Y.  Cheng,  and  K.  S.  Leung,  "Mechanical  properties  of  normal  skin  and 
hypertrophic  scars,"  Burns  22,  443-446  (1996). 

11.  H.  Dobrev,  "In  vivo  study  of  skin  mechanical  properties  in  psoriasis  vulgaris,"  Acta  Derm.- 
Venereol.  80,  263-266  (2000). 

12.  H.  Sumino,  S.  Ichikawa,  M.  Abe,  Y.  Endo,  O.  Ishikawa,  and  M.  Kurabayashi,  "Effects  of  aging, 
menopause,  and  hormone  replacement  therapy  on  forearm  skin  elasticity  in  women,"  J.  Am.  Geriatr.  Soc. 
52,945-949(2004). 

13.  A.  A.  Marghoob,  L.  D.  Swindle,  C.  Z.  M.  Moricz,  F.  A.  S.  Negron,  B.  Slue,  A.  C.  Halpem,  and 
A.  W.  Kopf,  "Instruments  and  new  technologies  for  the  in  vivo  diagnosis  of  melanoma,"  Journal  of  the 
American  Academy  of  Dermatology  49,  777-797  (2003). 

14.  E.  Ruocco,  G.  Argenziano,  G.  Pellacani,  and  S.  Seidenari,  "Noninvasive  imaging  of  skin  tumors," 
Dermatol.  Surg.  30,  301-310  (2004). 

15.  R.  Dua,  D.  G.  Beetner,  W.  V.  Stoecker,  and  D.  C.  Wunsch,  "Detection  of  basal  cell  carcinoma 
using  electrical  impedance  and  neural  networks,"  IEEE  Trans.  Biomed.  Eng.  51,  66-71  (2004). 

16.  M.  Gniadecka,  P.  A.  Philipsen,  S.  Sigurdsson,  S.  Wessel,  O.  F.  Nielsen,  D.  H.  Christensen,  J. 
Hercogova,  K.  Rossen,  H.  K.  Thomsen,  R.  Gniadecki,  L.  K.  Hansen,  and  H.  C.  Wulf,  "Melanoma 
diagnosis  by  Raman  spectroscopy  and  neural  networks:  Structure  alterations  in  proteins  and  lipids  in 
intact  cancer  tissue,"  J.  Invest.  Dermatol.  122,  443-449  (2004). 

17.  C.  Edwards,  and  R.  Marks,  "Evaluation  of  biomechanical  properties  of  human  skin,"  Clin. 
Dermatol.  13,  375-380  (1995). 

18.  M.  I.  Miga,  "A  new  approach  to  elastography  using  mutual  information  and  finite  elements," 
Phys.  Med.  Biol.  48,  467-480  (2003). 

19.  L.  V.  Tsap,  D.  B.  Goldgof,  and  S.  Sarkar,  "Nonrigid  motion  analysis  based  on  dynamic 
refinement  of  finite  element  models,"  IEEE  Trans.  Pattern  Anal.  Mach.  Intell.  22,  526-543  (2000). 

20.  L.  V.  Tsap,  D.  B.  Goldgof,  and  S.  Sarkar,  "Fusion  of  physically-based  registration  and 
deformation  modeling  for  nonrigid  motion  analysis,"  IEEE  Trans.  Image  Process.  10,  1659-1669  (2001). 

21.  L.  V.  Tsap,  D.  B.  Goldgof,  S.  Sarkar,  and  P.  S.  Powers,  "A  vision-based  technique  for  objective 
assessment  of  burn  scars,"  IEEE  Trans.  Med.  Imaging  17,  620-633  (1998). 

22.  T.  L.  Chenevert,  A.  R.  Skovoroda,  M.  O'Donnell,  and  S.  Y.  Emelianov,  "Elasticity  reconstructive 
imaging  by  means  of  stimulated  echo  MRI,"  Magn.Reson.Med.  39,  482-490  (1998). 

23.  J.  Ophir,  S.  K.  Alam,  B.  Garra,  F.  Kallel,  E.  Konofagou,  T.  Krouskop,  and  T.  Varghese, 
"Elastography:  ultrasonic  estimation  and  imaging  of  the  elastic  properties  of  tissues,"  Proc.  Inst.  Mech. 
Eng.  Part  H-J.  Eng.  Med.  213,  203-233  (1999). 


35 


24.  E.  E.  W.  Van  Houten,  M.  I.  Miga,  J.  B.  Weaver,  F.  E.  Kennedy,  and  K.  D.  Paulsen,  "Three- 
dimensional  subzone-based  reconstruction  algorithm  for  MR  elastography,"  Magn.Reson.Med.  45,  827- 
837  (2001). 

25.  M.  I.  Miga,  M.  P.  Rothney,  and  J.  J.  Ou,  "Modality  independent  elastography  (MIE):  Potential 
applications  in  dermoscopy,"  Medical  Physics  32,  1308-1320  (2005). 

26.  M.  I.  Miga,  "A  new  approach  to  elastographic  imaging:  Modality  independent  elastography," 
Medical  Imaging  2002:  Image  Processing:  Proc.  of  the  SPIE  4684,  604-61 1  (2002). 

27.  J.  J.  Ou,  S.  L.  Barnes,  and  M.  I.  Miga,  "Application  of  multi-resolution  modality  independent 
elastography  for  detection  of  multiple  anomalous  objects,"  in  Medical  Imaging  2006:  Physiology, 
Function,  and  Structure  from  Medical  Images:  Proc.  of  the  SPIE,  (SPIE,  San  Diego,  CA,  2006). 

28.  J.  J.  Ou,  S.  L.  Bames,  and  M.  I.  Miga,  "Preliminary  testing  of  sensitivity  to  input  data  quality  in 
an  elastographic  reconstruction  method,"  in  International  Symposium  on  Biomedical  Imaging  2006, 
(Washington,  D.C.,  2006). 

29.  C.  W.  Washington,  and  M.  I.  Miga,  "Modality  independent  elastograph  ({MIE}):  {A}  new 
approach  to  elasticity  imaging,"  IEEE  Trans.  Med.  Imaging  (in  press),  (2004). 

30.  D.  W.  Marquardt,  "An  algorithm  for  least-squares  estimation  of  nonlinear  parameters,"  SIAM 
Journal  of  Applied  Mathematics  11,  431-441  (1963). 

31.  N.  Joachimowicz,  C.  Pichot,  and  J.  P.  Hugonin,  "Inverse  scattering:  an  iterative  numerical 
method  for  electromagnetic  imaging,"  IEEE  Trans.  Biomed.  Eng.  39,  1742-1752  (1991). 

32.  R.  Horejsi,  R.  Moller,  T.  R.  Pieber,  S.  Wallner,  K.  Sudi,  G.  Reibnegger,  and  E.  Tafeit, 
"Differences  of  subcutaneous  adipose  tissue  topography  between  type-2  diabetic  men  and  healthy 
controls,"  Exp.  Biol.  Med.  227,  794-798  (2002). 

33.  R.  Moller,  E.  Tafeit,  K.  H.  Smolle,  T.  R.  Pieber,  O.  Ipsiroglu,  D.  M.,  C.  Huemer,  K.  Sudi,  and  G. 
Reibnegger,  "Estimating  percentage  total  body  fat  and  determining  subcutaneous  adipose  tissue 
distribution  with  a  new  non-invasive  optical  device  LIPOMETER,"  Americal  Journal  of  Human  Biology 
12,221-230(2000). 

34.  F.  M.  Hendriks,  D.  Brokken,  J.  T.  W.  M.  van  Eemeren,  C.  W.  J.  Oomens,  F.  P.  T.  Baaijens,  and  J. 
B.  A.  M.  Horsten,  "A  numerical-experimental  method  to  characterize  the  non-linear  mechanical  behavior 
of  human  skin,"  Skin  Res.  Technol.  9,  274-283  (2003). 


36 


Skin  Cancer  &  Significance 


An  elastography  framework 
for  use  in  dermoscopy 


Michael  I.  Miga,  Jao  J.  Ou,  Darrel  L.  Ellis 
Vanderbilt  University 
Vanderbilt  University  Medical  Center 
February  20,  2007 

V _ _ _ / 

o 

bml 


•  Typical  Discerned  Facts  for  Skin  Cancer: 

3  types  -  basal  cell  carcinoma  (BCC),  squamous  cell 

carcinoma  (SCC),  and  melanoma 

1,000,000+  new  cases  are  diagnosed  each  year 

800,000+  BCCs,  200,000+  SCCs,  60,000+  melanomas 

Annually,  melanoma  ranks  6th  highest  among  newly 

diagnosed  cancers 

Melanoma  is  most  lethal 

Skin  lesion  scoring  for  monitoring  (ABCDE,  etc.) 

Diagnosis  is  performed  by  biopsy 

Skin  cancer  is  rarely  fatal  and  is  largely  treatable 


\ 


If  so,  why  pursue  this? 


Skin  Cancer  &  Significance 
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•  Less  Discerned  Facts  for  Skin  Cancer: 

1,000,000+  newly  diagnosed  skin  cancers  involve  an 
estimated  5-7  million  biopsies 
Management  of  NMSC  is  considerable 


Housman  et  al.,  Journal  of  the  American  Academy  of  Dermatology, 
vol.  48,  pp.  425-429,  Mar  2003. 


•  Green  et  al.,  Dermatologic  Surgery,  vol.  30,  pp.  1208- 
1209,  Sep  2004. 

22  Dermatologists  participated  to  determine  how 
many  biopsies  taken  were  cancerous 
“Malignancy  ratio”  (#  of  malignant  biopsies  to  total 
biopsies  performed) 

Findings: 

•  40%  mean  malignancy  ratio  (ratios  varied  greatly  13-85%) 
Age  correlated  with  better  malignancy  ratio’s  (2%  increase 
with  year  of  experience) 


Annual  biopsy  costs  are  in  the  $1 00 ’s  of  millions 


Motivation 


•  Total  direct  cancer  costs  are  approximately  78.2  billion 
per  year 

•  If  a  low-cost  dermoscopic  probe,  with  elastographic 
imaging  capabilities  was  to  have  high  sensitivity  and 
specificity  for  skin  cancers,  this  would  translate  to  a 
significant  reduction  in  health  care  expenditure. 


•  While  skin  cancer  may  lack  significance  in 
lethality,  it  is  quite  significant  with  respect  to 
healthcare  costs. 


Elasticity  Imaging  Approach 


•  Computational  Modeling 

Patient-specific  finite  element  models 
Physics-based  simulations  and  image  deformation 

•  Image  Similarity 

Metrics  of  image  comparison 

•  Non-linear  Optimization 

Cornerstone  in  model-based  image  reconstruction 

methods 

Inverse  problem 

A  Modality  Independent  Elastography(MIE)J 
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MIE  Method  for  Dermoscopic  Data  Sets 


Paper  Purpose 


•  Investigate  the  use  of  a  geometric 
nonlinear  model  for  inversion 


•  Investigate  the  effects  of  3D  nature  of 
skin  (both  lesion  and  subcutaneous 
layers) 


Present  a  very  preliminary  in  vivo  result 


Elasticity  Image  Reconstructions 


Bulk 

Experiment 

Experiment 

Mean  Reconstructed 

Mean  Reconstructed 

Material 

Contrast  Ratio 

Contrast  Ratio 

Contrast  Ratio 

Contrast  Ratio 

Stretch 

Linear  Strain 

NonLinear  Strain 

Linear  Model  (Max) 

NonLinear  Model  (Max) 

5  mm 

5.7 

5.0 

3.8+/-0.8  (5.0) 

3.5+/-0.6  (4.4) 

10  mm 

5.0 

4.8 

3.2+/-0.7  (4.2) 

3.5+/-0.8  (4.5) 

15  mm 

4.6 

4.7 

3.5+/-0.9  (5.1) 

3.7+/-0.9  (5.1) 

20  mm 

4.2 

4.4 

2.8+/-0.7  (3.8) 

2.9+/-0.8  (4.2) 

Elasticity  Image  Reconstructions 
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Paper  Purpose 


Simulation  Setup 
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•  Investigate  the  use  of  a  geometric 
nonlinear  model  for  inversion 


•  Investigate  the  effects  of  3D  nature  of 
skin  (both  lesion  and  subcutaneous 
layers) 


Present  a  very  preliminary  in  vivo  result 


•  Systematic  stretch 

•  Dermoscopic  video 
image  acquisition 


Results 
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Paper  Purpose 


•  Investigate  the  use  of  a  geometric 
nonlinear  model  for  inversion 


•  Investigate  the  effects  of  3D  nature  of 
skin  (both  lesion  and  subcutaneous 
layers) 

•  Present  a  very  preliminary  in  vivo  result 


Investigating  In  Vivo 


•  Spring-loaded  grips  bounded  to  the  skin 
with  adhesive 
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Other  Membrane  Exoeriments 


Conclusions 


Geometric  nonlinear  effects  are  mild 

Boundary  conditions  are  likely  to  be  the  most 
important  aspect 

5-10%  improvement  in  objective  function  value 
Qualitative  improvement  in  lesion  border 
Material  nonlinearities  are  next  step 
Subcutaneous  layers  have  little  effect 
Depth  differentiation  with  contrast 
Overlap  in  contrast  may  be  problem 
In  vivo  results  indicate  a  nevus  lesion  has 
increased  stiffness 


Future  Work 


•  Device 

•  MIE’s  robustness 

•  Integration 

•  3D  Extension 
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to  Boundary  Condition  Noise 
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ABSTRACT 

This  work  explores  an  inverse  problem  technique  of  extracting  soft  tissue  elasticity  information  via  nonrigid  model- 
based  image  registration.  The  algorithm  uses  the  elastic  properties  of  the  tissue  in  a  biomechanical  model  to  achieve 
maximal  similarity  between  image  data  acquired  under  different  states  of  loading.  A  framework  capable  of  handling 
fully  three-dimensional  models  and  image  data  has  been  recently  developed  utilizing  parallel  computing  and  iterative 
sparse  matrix  solvers.  For  this  preliminary  investigation,  a  series  of  simulation  experiments  with  clinical  image  data  of 
human  breast  are  used  to  test  the  robustness  of  the  algorithm  to  expected  mis-estimation  of  displacement  boundary 
conditions  encountered  in  real-world  situations.  Three  methods  of  automated  point  correspondence  are  also  examined  as 
means  of  generating  boundary  conditions  for  the  algorithm. 

Keywords:  elastography,  computational  modeling,  inverse  problem,  non-rigid  registration 


1.  INTRODUCTION 

The  characterization  of  the  mechanical  properties  of  tissue  is  an  important  potential  source  of  information  for  detection 
and  diagnosis  of  disease  processes.  For  example,  there  is  a  long-standing  clinical  appreciation  of  evaluating  tissue 
elasticity  through  palpation  in  the  physical  examination  and  correlating  differences  in  stiffness  with  possible 
pathological  states.  A  minimally  invasive  methodology  for  analyzing  tissue  deformation  through  imaging  and/or  image 
processing  techniques  is  a  central  goal  of  the  field  of  elastography  [1,11].  Application  of  such  methods  to  the 
interrogation  of  the  breast  [2,3],  skin  [4-6],  prostate  [7],  and  other  accessible  organ  systems  is  an  active  area  of  research. 

Many  of  the  current  elastography  methods  are  founded  in  ultrasound  (US)  and  magnetic  resonance  (MR)  imaging  and 
involve  the  estimation  of  induced  displacements  within  the  tissue  of  interest  to  infer  the  elasticity  distribution.  We  have 
recast  the  problem  as  a  physically-constrained  non-rigid  image  registration  utilizing  quasi-static  deformation  and  image 
similarity  metrics  that  reconstruct  the  spatial  distribution  of  elasticity  parameters.  This  technique  has  been  termed 
’modality-independent  elastography’  (MIE)  [8-10]  because  of  its  ability  to  handle  native  anatomical  images  from 
different  sources  with  relatively  simple  modifications  to  the  acquisition  procedure.  To  date,  data  from  MR,  X-ray 
computed  tomography  (CT),  and  digital  photography  have  been  used  to  drive  the  algorithm.  In  addition  to  the  necessary 
preparation  and  effort  involved  in  gathering  images,  the  other  major  input  to  this  reconstruction  method  is  the 
delineation  of  boundary  conditions  on  the  region  of  interest.  Because  this  process  currently  involves  varying  levels  of 
manual  interaction,  there  is  a  need  to  develop  a  protocol  that  is  both  effective  and  mostly  automated  for  determining 
point  correspondences.  The  objectives  of  this  work  are  to  test  the  effects  of  degradation  in  input  data  quality  on  the  end 
reconstruction  and  examine  candidate  methods  for  generation  of  displacement  boundary  conditions.  This  is  done  in  the 
context  of  evaluating  the  robustness  of  a  newly  realized  three-dimensional  version  of  MIE  by  performing  simulation 
experiments  with  randomized  noise  processes  and  comparing  the  fidelity  of  reconstructions  resulting  from  boundary 
conditions  generated  by  three  different  techniques  of  determining  surface  point  correspondence. 

2.  METHODS 


2.1  Elastographic  reconstruction  framework 

The  conceptual  framework  for  our  elastographic  reconstruction  has  been  previously  described  in  [6,8-10].  In  brief,  an 
image  of  a  tissue  of  interest  {source)  is  deformed  by  a  biomechanical  computer  model  and  compared  against  an  acquired 
image  of  the  same  tissue  in  a  mechanically  loaded  state  {target).  Iterative  updates  of  elasticity  parameters  to  the  model 
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are  performed  until  a  suitable  match  in  intramodal  image  similarity  is  achieved  in  a  least  squares  manner  to  satisfy  a 
non-linear  optimization  scheme.  This  process  can  be  classified  as  an  inverse  problem,  with  model-based  deformation  of 
the  source  image  representing  the  forward  problem.  The  three  major  components  of  the  algorithm  are  the  model,  image 
comparison,  and  optimization,  each  of  which  is  described  in  more  detail  below. 

The  partial  differential  equation  that  expresses  a  state  of  mechanical  equilibrium  can  be  written  as  [12]: 

V  -cr  =  0  (1) 


where  <ris  the  Cartesian  stress  tensor.  We  have  elected  to  model  the  constitutive  tissue  behavior  using  Hooke’s  Law  of 
linear  elasticity,  which  states  that  the  strain  is  proportional  to  the  applied  stress,  and  assume  that  materials  are  isotropic 
and  incompressible  in  nature.  This  leads  to  the  formulation  of  Cauchy’s  Law  as 
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which  describes  the  constitutive  relationship  between  stress  and  strain  in  terms  of  the  elasticity  parameters  E  (Young’s 

£ 

modulus)  and  v  (Poisson’s  ratio).  The  shear  modulus  G  is  defined  as - . 

2(1  + v) 


A  finite  element  (FE)  representation  of  the  model  is  constructed  from  the  source  image  and  assigned  appropriate 
boundary  conditions  based  on  estimated  displacement  or  stress.  The  standard  Galerkin  method  of  weighted  residuals 
[13]  is  used  to  construct  and  solve  the  system,  which  yields  a  set  of  displacements  that  are  used  to  deform  the  source 
image.  This  model-deformed  image  is  then  compared  to  the  target  using  an  intensity-based  image  similarity  calculated 
for  a  series  of  voxel  groupings  determined  by  a  downsampling  of  the  image  set  overlap.  The  correlation  coefficient 
(CC)  [14]  is  the  method  of  choice,  as  it  has  empirically  demonstrated  superior  performance  over  other  metrics  such  as 
the  sum  of  squared  differences  and  normalized  mutual  information. 

The  elasticity  parameter  optimization  can  be  written  as  the  minimization  of  the  least  squares  error  objective  function 


Y  =  \s. 


TRUE 


f  EST 


(3) 


where  STrue  is  the  set  of  similarity  values  achieved  when  comparing  the  target  image  to  itself,  SEst  is  the  similarity 
between  the  model-deformed  source  and  the  target  images  using  current  estimates  of  the  elastic  modulus  distribution, 
and  |*|  denotes  the  vector  L2  norm.  Note  that  by  definition,  STRUE  for  CC  has  a  constant  value  of  1.  Using  a  Levenberg- 
Marquardt  approach,  the  residual  form  of  equation  (3)  becomes 


[jr/  +  a/]!A£}=[jr|sr„l,I  -Slsr} 
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where  J  =  5SEst/^E  is  the  Jacobian  matrix  and  the  regularization  parameter  a  is  determined  using  the  methods  described 
in  [15].  Modulus  values  are  updated  by  AE  until  an  error  tolerance  is  reached  or  a  maximum  number  of  iterations  have 
been  completed.  Spatial  averaging  of  elasticity  values  in  the  model  and  solution  relaxation  between  iterations  are  also 
utilized  to  improve  the  stability  of  the  optimization. 

It  should  be  noted  that  the  size  of  the  Jacobian  matrix  is  dependent  on  the  number  of  material  properties  to  be 
reconstructed,  with  each  column  requiring  a  forward  solve  of  the  FE  model.  For  the  general  lesion  detection  problem,  a 
fine  discretization  of  the  mesh  requires  many  solutions  such  that  the  building  of  this  matrix  consumes  a  considerable 
portion  of  computational  resources.  This  fact  is  exacerbated  with  the  use  of  three-dimensional  data  and  necessitates  a 
parallelized  system.  Recent  work  using  the  Portable  Extensible  Toolkit  for  Scientific  Computation  (PETSc)  toolkit 
[16,17]  has  provided  the  necessary  coding  base  for  interfacing  sparse  matrix  system  solvers  with  a  C/C++  optimization 
framework  in  order  to  supersede  the  original  MATLAB/FORTRAN  framework.  This  new  version  is  designed  to  scale 
readily  between  the  complexity  of  the  model  and  the  number  of  processors  available;  it  has  been  tested  on  a 
homogeneous  cluster  of  ten  processors,  with  further  active  development  taking  place  in  conjunction  with  the  Vanderbilt 
ACCRE  project,  which  houses  hundreds  of  CPUs. 

2.2  Simulation  experiment  setup 

A  CT  volume  of  a  human  breast,  obtained  from  UC-Davis  Dept,  of  Radiology,  was  used  as  the  source  image  (256  x  256 
x  130,  0.6mm  x  0.6mm  x  0.6mm  voxel  spacing)  for  the  remainder  of  this  work.  The  surface  of  the  breast  was  segmented 
(ANALYZE  6.0,  Mayo  Cline,  Rochester,  MN)  to  create  a  three-dimensional  mesh  composed  of  39,013  nodes  connected 
as  214,163  tetrahedral  elements.  In  order  to  ensure  initial  data  fidelity  for  reconstruction  experiments,  an  idealized  target 
image  volume  and  gold  standard  boundary  condition  set  were  created.  A  2-cm  spherical  tumor  was  implanted  in  the 
center  of  the  mesh  by  assigning  a  stiff  modulus  value  to  the  member  elements  that  was  six  times  higher  than  the 
surrounding  material  [18].  Tissue  deformation  from  the  inflation  of  a  rectangular  air  bladder  against  the  lateral  surface  of 
the  breast  was  numerically  simulated  to  qualitatively  match  observed  mechanically  loading  of  an  existing  device  on  a 
breast-mimicking  phantom  of  polyvinyl  alcohol  cryogel.  The  stress  distribution  over  a  rectangular  contact  area  was 
modeled  as  the  cross-section  of  a  Gaussian  pressure  field  with  its  maximum  value  located  at  the  center  of  the  bladder; 
the  base  of  the  breast  was  fixed  in  place  as  if  pinned  to  the  chest  wall.  The  deformation  field  throughout  the  domain  was 
calculated  using  a  direct  forward  solve  of  the  finite  element  model  and  then  applied  to  the  intensity  field  of  the  source 
image  to  create  a  target  volume.  Displacements  at  the  surface  nodes  were  used  to  make  a  final  description  comprised  of 
all  Type  I  (Dirichlet)  boundary  conditions.  Figure  1  below  illustrates  the  setup  of  the  data  used. 

All  reconstructions  were  performed  using  a  priori  knowledge  of  the  location  and  size  of  the  inclusion  in  order  to  limit 
the  scope  of  the  problem  (e.g.  the  expense  of  the  Jacobian  matrix)  to  a  two-material  discrimination  of  relative  stiffness 
(elastic  contrast).  Having  a  defined  physical  model  and  synthetic  image  comparison  also  allows  for  examination  of  the 
optimization  behavior  separately  from  the  other  MIE  components  in  order  to  best  evaluate  the  effect  of  input 
inaccuracies  on  the  final  elasticity  distribution.  The  reconstruction  algorithm  begins  by  assigning  a  homogeneous 
elasticity  distribution,  with  Poisson’s  ratio  held  constant  at  v  =  0.485  to  represent  a  nearly  incompressible  material.  For 
this  data  set,  733  similarity  zones  were  discretized  from  the  target  image  volume. 


Figure  1.  Surface  renderings  of  CT  breast  volume  used  in  MIE  simulation  experiments.  From  left  to  right:  source  image, 
finite  element  mesh,  and  target  image  with  deformation  created  by  presumed  inflation  of  an  air  bladder. 
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2.3  Testing  robustness  of  the  algorithm 

The  current  method  of  selecting  boundary  conditions  as  derived  from  experiences  in  2D  work  requires  manual 
interaction  to  guide  or  correct  point  correspondence  for  every  surface  node.  Assuming  that  visible  markers  are  available 
in  an  image,  but  that  an  input  device  (e.g.  a  mouse)  is  needed  to  identify  the  specific  coordinates,  this  introduces  an 
operator-dependent  noise  process  in  localizing  any  given  point.  The  cumulative  effect  of  these  inaccuracies  is 
propagated  from  the  model  to  the  image  deformation  and  then  the  similarity  measurements.  For  a  given  reconstruction 
experiment,  the  gold  standard  boundary  condition  set  was  systematically  disrupted  by  adding  a  Gaussian  randomized 
vector  of  a  particular  length  (0.1,  0.2,  0.5,  1.0  or  2.0  voxel  units).  Figure  2  shows  an  example  of  the  distortion  caused  by 
the  applied  noise. 


Figure  2.  Example  of  distortion  due  to  additive  randomized  error.  The  gold  standard  boundary  conditions  used  to  generate 
a  controlled  deformation  produce  the  mesh  shown  on  the  left.  For  effect,  the  2.0  voxel  unit  noise  is  show  on  the  right. 


2.4  Testing  automated  boundary  condition  generating  methods 

For  this  work,  three  methods  of  surface  registration  and  point  correspondence  were  considered  for  a  more  automated 
method  of  determining  boundary  conditions  for  the  reconstruction  algorithm.  Two  are  derived  from  surface  matching  of 
potential  energy  distributions,  and  the  other  is  a  free-form  warping. 

If  the  flow  of  a  substance  over  undeformed  and  deformed  breast  surfaces  is  taken  to  be  a  conserved  process,  then 
correspondence  can  be  achieved  by  matching  the  same  energy  deposition  between  the  source  and  target,  that  is,  the 
equivalent  level  sets.  In  a  physical  sense,  Laplace’s  equation  can  be  used  to  describe  this  type  of  movement  analogously 
to  steady-state  heat  distribution: 

V2®  =  0  (5) 


where  O  would  represent  the  temperature  over  a  given  region.  Similarly,  the  diffusion  equation  describes  the  change  in 
concentration  or  density  of  a  material  over  time  on  a  region: 

5®  „  o -r. 

—  =  aV  ®  (6) 

dt 


where  a  is  the  diffusivity  constant. 

These  partial  differential  equations  were  used  to  calculate  an  energy  distribution  from  the  nipple  area  to  the  chest  wall 
over  the  surface  of  a  breast.  Isocontours  of  particular  energy  values  were  then  extracted  from  each  surface  to  form  a  set 
of  connected  points.  The  symmetric  closest  point  method  described  by  [19]  was  used  to  determine  a  displacement  field 
from  which  point  correspondence  at  boundary  nodes  could  be  interpolated. 
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The  third  method  involves  thin-plate  spline  interpolation  [20]  to  determine  point  correspondence.  This  was  done  to 
consider  a  widely-used  method  of  non-rigid  transform  that  can  take  advantage  of  fiducial  information  that  should  be 
present  in  future  real-world  data  acquisition.  The  use  of  physical  markers  to  track  breast  surface  displacement  during  a 
deformation  also  defines  a  set  of  control  points  would  allow  the  displacement  boundary  conditions  for  MIE  to  be  simply 
interpolated  from  the  local  warping.  In  these  simulation  experiments,  40  surface  nodes  of  known  correspondence  in 
each  of  the  image  volumes  (due  to  the  controlled  deformation)  served  as  fiducials. 

2.5  Evaluation  of  reconstructions 

Evaluation  of  the  reconstruction  results  is  performed  by  calculating  the  ratio  of  the  elasticity  of  the  inclusion  to  the  rest 
of  the  breast  for  the  distribution  that  yields  the  minimal  objective  function  value  over  the  course  of  optimization.  The 
robustness  of  the  MIE  algorithm  was  tested  with  four  trials  at  each  of  the  magnitudes  of  randomized  vectors  (described 
above  in  Section  2.3),  and  the  reconstruction  results  were  averaged  to  determine  a  trend  and  possible  threshold  of  noise 
tolerance  for  the  algorithm. 

For  the  automated  boundary  condition  generation  methods,  a  forward  mapping  of  the  objective  function  space  was 
calculated  to  determine  a  theoretical  optimum  to  the  reconstruction.  This  was  done  by  calculating  the  similarity  values 
for  model-based  image  deformations  created  by  adjusting  the  elasticity  contrast  of  the  inclusion  over  a  range  of  0.5:1  to 
30:1.  An  interpolating  curve  was  fit  and  the  minimum  objective  function  value  and  associated  elasticity  contrast  were 
extracted. 


3.  RESULTS 


3.1  Robustness  of  algorithm  to  boundary  condition  noise 

The  following  tables  summarize  the  effects  of  additive  noise  of  a  particular  magnitude  to  the  gold  standard  boundary 
condition  set.  As  the  magnitude  of  the  applied  randomized  vectors  increased,  a  dramatic  increase  in  the  minimum 
objective  function  value  is  observed.  Additionally,  changes  in  the  reconstructed  elasticity  contrast  indicate  that  a  cutoff 
exists  in  the  ability  of  the  algorithm  to  achieve  a  successful  result  (recall  that  the  known  correct  ratio  is  6:1)  for 
disruption  by  vectors  of  length  0.5  voxel  units  or  higher. 


Table  1.  Effect  of  applied  random  boundary  condition  noise  on  objective  function  space  and  reconstructed  elasticity 
contrast  ratio. 


Randomized 
vector  magnitude 
(voxel  units) 

Mean  optimal 
objective  function 
value 

Mean  optimal 
elasticity  contrast 
value 

0.1 

2.85  ±0.0382 

5.62  ±0.421 

0.2 

10.1  ±0.367 

5.70  ±0.588 

0.5 

60.1  ±4.19 

2.36  ±0.393 

1.0 

80.2  ±0.561 

2.47  ±  0.266 

2.0 

104  ±3.42 

2.17  ±0.422 
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3.2  Use  of  automated  point  correspondence 

Figure  3  depicts  the  deformation  fields  found  by  the  various  automated  methods  which  were  converted  into  a  boundary 
condition  sets  and  run  through  the  reconstruction  algorithm.  Qualitatively,  the  displacements  found  by  the  diffusion 
method  are  quite  different  from  the  gold  standard  set,  while  the  results  from  the  solution  of  Laplace’s  equation  and  the 
thin-plate  spline  interpolation  appear  to  be  more  satisfactory.  The  mean  target  registration  error  of  the  three  methods 
confirms  this  with  the  spline  having  the  best  performance  (0.26  mm),  the  Laplace  method  next  (1.0  mm),  and  the 
diffusion  method  being  the  worst  (2.0  mm).  Inspection  of  Figure  4  further  demonstrates  that  the  imposition  of  an  inexact 
boundary  condition  set  on  the  model  has  a  distinct  effect  on  the  optimization  by  shifting  the  minimum  value  to  a 
different  position.  A  comparison  of  the  fit  with  the  reconstruction  in  both  objective  function  value  and  elasticity  contrast 
is  provided  in  Table  3  below  and  indicates  that  the  algorithm  is  mostly  in  agreement  with  the  predicted  values  for  the 
Laplace  and  thin-plate  spline  methods  but  not  as  well  for  the  diffusion  method. 


Figure  3.  Three  candidate  automated  methods  for  MIE  boundary  condition  generation.  Top  row,  from  left  to  right:  surface 
deformations  calculated  from  diffusion  energy  matching,  Laplace  solution  energy,  and  thin-plate  spline  interpolation. 
Bottom  row:  target  registration  error  (TRE)  distribution  for  each  method  when  compared  against  the  gold  standard  of 
known  correspondence.  The  diffusion-based  mesh  is  both  qualitatively  and  quantitatively  the  worst  performer.  The 
Laplace  solution  appears  to  capture  the  shape  of  the  bladder  indentation  more  precisely,  but  the  thin-plate  spline  has  the 
best  overall  accuracy  in  determining  point  correspondence. 
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Figure  4.  Mappings  of  objective  function  value  vs.  elasticity  contrast  ratio  (tumoribreast)  as  affected  by  the  boundary 
condition  sets  generated  by  different  automated  methods  of  surface  point  correspondence.  The  minimum  value  of  each 
curve  corresponds  to  the  optimal  elasticity  contrast  that  can  be  achieved  by  the  algorithm  when  constrained  by  the 
inaccuracies  of  the  methods:  (a)  diffusion,  (b)  Laplace,  and  (c)  thin-plate  spline  interpolations. 


Table  3.  Comparison  of  automated  point  correspondence  methods  on  MIE  reconstruction  quality.  Predicted  values  are 
found  from  the  minimum  point  of  the  curves  shown  in  Figure  4. 


Method 

Predicted  minimum 
objective  function 
value 

Reconstructed 
minimum  objective 
function  value 

Predicted 
optimal  elasticity 
contrast 

Reconstructed 
optimal  elasticity 
contrast 

Diffusion 

58.5 

58.8 

30.0 

12.6 

Laplace 

5.59 

5.59 

9.55 

10.3 

Thin-plate  spline 

12.3 

12.3 

5.55 

5.66 
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4.  DISCUSSION 


The  results  of  the  boundary  condition  noise  experiment  are  interesting  because  they  indicate  that  improper  localization 
of  boundary  points  greater  than  or  equal  to  0.5  units  of  voxel  spacing  can  introduce  significant  error  to  the  reconstruction 
process  and  impair  its  ability  to  characterize  the  underlying  elasticity  distribution.  This  is  a  similar  result  to  prior  work 
done  in  two-dimensional  systems  in  which  successful  reconstructions  correlated  to  boundary  condition  selection  error 
limited  to  half  a  pixel  length  [21].  It  also  confirms  that  randomizing  the  vectors  is  a  significant  challenge  to  the  algorithm 
because  it  introduces  highly  non-physical  deformations  that  cause  backlash  in  the  finite  element  mesh  and  other 
numerical  anomalies. 

The  implausibility  of  performing  manual  selection  on  all  6,319  boundary  nodes  underscores  the  importance  of  finding  an 
automated  method  for  determining  point  correspondences,  especially  at  less  than  0.5  voxel  units  of  error.  In  these 
simulation  experiments,  energy  matching  from  the  solutions  of  the  diffusion  and  Laplace  equations  yield  boundary 
condition  sets  that  are  inadequate  for  reconstructing  a  proper  elasticity  contrast.  This  can  be  partly  explained  because  the 
mean  errors  of  those  surface  registration  techniques  (as  compared  to  the  gold  standard)  are  approximately  3.3  and  1.7 
voxel  units,  respectively,  which  based  on  the  randomized  trials  were  magnitudes  too  large  for  the  algorithm  to  handle. 
The  diffusion-based  boundary  conditions  also  proved  more  difficult  to  obtain  a  stable  solution  for  in  the  model,  which 
probably  contributed  to  the  mismatch  in  reconstructed  elasticity  contrast.  However,  the  results  obtained  from 
reconstructions  using  the  thin-plate  spline  method  are  encouraging  because  the  mean  error  was  0.43  voxel  units.  The 
reconstruction  behavior  in  that  case  was  consistent  with  the  predicted  objective  function  space  and  the  optimal  elasticity 
contrast  was  found  to  be  within  6%  of  the  true  value.  This  preliminary  result  appears  to  identify  the  use  of  thin-plate 
spline  interpolation  as  a  strong  candidate  for  generating  boundary  conditions  for  MIE.  The  use  of  40  control  points  is 
also  seen  as  a  reasonable  choice  for  data  acquisition  and  processing  in  order  to  capture  the  extent  of  expected 
deformation  processes. 


5.  CONCLUSIONS 

In  this  work,  we  have  demonstrated  the  effects  of  inaccuracies  in  boundary  condition  determination  on  an  elastography 
method  that  maximizes  the  similarity  between  images  of  a  tissue  of  interest  acquired  under  two  different  states  of 
mechanical  loading.  In  order  to  characterize  the  robustness  of  the  current  version  of  this  method,  which  has  been 
updated  to  handle  three-dimensional  data  in  a  parallelized  scheme,  randomized  vectors  were  applied  to  distort  a  gold 
standard  boundary  condition  set.  The  results  were  used  to  determine  a  threshold  of  accuracy  needed  in  order  to  still 
achieve  an  accurate  reconstruction.  In  order  to  streamline  the  pre-processing  involved  in  the  algorithm,  three  methods  of 
automated  point  correspondence  were  evaluated.  The  success  of  these  methods  correlated  with  their  mean  error  (relative 
to  the  true  displacements)  meeting  the  putative  cutoff,  and  initial  results  indicate  that  established  techniques  such  as  thin- 
plate  splines  hold  promise  for  generating  boundary  conditions. 
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BACKGROUND 


This  work  explores  an  inverse  problem  technique  of  extracting  soft  tissue  elasticity  information  via  nonrigid  model-based 
image  registration.  The  algorithm  uses  the  elastic  properties  of  the  tissue  in  a  biomechanical  model  to  achieve  maximal 
similarity  between  image  data  acquired  under  different  states  of  loading.  A  framework  capable  of  handling  fully  three- 
dimensional  models  and  image  data  has  been  recently  developed  utilizing  parallel  computing  and  iterative  sparse  matrix 
solvers.  For  this  preliminary  investigation,  a  series  of  simulation  experiments  with  clinical  image  data  of  human  breast  are 
used  to  test  the  robustness  of  the  algorithm  to  expected  mis-estimation  of  displacement  boundary  conditions  encountered  in 
real-world  situations.  Three  methods  of  automated  point  correspondence  are  also  examined  as  means  of  generating 
boundary  conditions  for  the  algorithm. 


METHODS 


The  conceptual  framework  for  the  elastographic  reconstruction  involves  three  major  components:  a  biomechanical  model  of 
the  tissue,  image  deformation  and  comparison,  and  numerical  optimization  of  image  similarity.  An  image  of  a  tissue  of 
interest  (source)  is  deformed  by  a  prescribed  computational  model  and  compared  against  an  acquired  image  of  the  same 
tissue  in  a  mechanically  loaded  state  (target).  The  deformation  and  comparison  is  repeated  using  systematic  updates  of 
elasticity  parameters  until  a  suitable  match  in  image  similarity  is  achieved  to  satisfy  a  non-linear  least  squares  objective 
function  [1],  Because  the  goal  is  to  determine  a  spatial  mapping  of  tissue  elasticity,  this  process  can  be  classified  as  an 
inverse  problem,  with  model-based  deformation  of  the  source  image  representing  the  forward  problem.  For  the  current 
version  of  MIE,  linear  elasticity  is  used  as  the  biomechanical  model,  the  correlation  coefficient  as  the  image  similarity 
metric,  and  a  regularized  Levenburg-Marquardt  algorithm  as  the  optimization  method.  The  constitutive  stress-strain 
relationship  of  the  model  is  solved  using  finite  elements  with  the  Galerkin  method  of  weighted  residuals,  and  Dirichlet 
boundary  conditions  are  selected  from  corresponding  points  in  the  image  pair. 


Figure  1.  MIE  reconstruction  algorithm  flowchart.  After  acquisition,  source  and  target  images  (A)  are  discretized  into  regions  and  zones, 
respectively.  The  reconstruction  process  involves  the  use  of  iteratively  updated  elastic  modulus  values  (B,  E)  to  drive  a  finite  element  model- 
based  image  deformation  (C)  until  the  best  match  is  found  (D). 


A  CT  volume  of  a  human  breast  (obtained  from  the  UC-Davis  Dept,  of  Radiology)  was  used  as  the  source  image  for  this 
work.  The  surface  of  the  breast  was  segmented  to  create  a  mesh  composed  of  39,013  nodes  connected  as  214,163 
tetrahedral  elements.  An  idealized  target  image  volume  and  gold  standard  boundary  condition  set  were  created  by 
implanting  a  2-cm  spherical  tumor  in  the  center  of  the  mesh  and  assigning  a  stiff  modulus  value  to  its  member  elements  that 
was  six  times  higher  than  the  surrounding  material.  Tissue  deformation  from  the  inflation  of  a  rectangular  air  bladder  against 
the  lateral  surface  of  the  breast  was  numerically  simulated  to  qualitatively  match  observed  mechanically  loading  behavior  of 
an  existing  device. 

The  current  method  of  selecting  displacement  (Type  1)  boundary  conditions  requires  manual  interaction  to  guide  or  correct 
point  correspondence  for  every  surface  node.  Assuming  that  an  input  device  (e.g.  a  mouse)  is  needed  to  identify  the  specific 
coordinates,  this  introduces  an  operator-dependent  noise  process  in  localizing  any  given  point.  The  cumulative  effect  of 
these  inaccuracies  is  propagated  from  the  model  to  the  image  deformation  and  then  the  similarity  measurements.  For  a 
given  reconstruction  experiment,  the  gold  standard  boundary  condition  set  was  systematically  disrupted  by  adding  a 
Gaussian  randomized  vector  of  a  particular  length  (0.1,  0.2,  0.5,  1.0  or  2.0  voxel  units).  Figure  2  shows  an  example  of  the 
distortion  caused  by  the  applied  noise,  and  Table  1  summarizes  their  effects  on  reconstructions  performed  with  a  priori 
knowledge  of  the  inclusion  location. 

The  implausibility  of  performing  manual  selection  on  all  6,319  surface  nodes  underscores  the  importance  of  developing  an 
automated  method  for  determining  boundary  conditions.  Three  methods  of  surface  registration  and  point  correspondence 
were  considered;  two  are  derived  from  surface  matching  of  potential  energy  distributions  based  on  the  diffusion  and  Laplace 
equations,  and  the  other  is  a  free-form  warping  via  a  thin-plate  spline. 

Table  1.  Effect  of  applied  random  boundary  condition  noise  on  objective 
function  space  and  reconstructed  elasticity  contrast  ratio  (ideal  6:1). 


Figure  2.  Top  row:  Surface  rendering  of  source  CT  breast  volume  (left)  and  finite  element  mesh  constructed  (right).  Bottom  row:  Target  image 
created  by  simulated  inflation  of  an  air  bladder  (left)  and  example  of  randomized  distortion  of  the  true  boundary  conditions. 


Randomized 
vector 
magnitude 
(voxel  units) 

Mean  optimal 
objective 
function  value 

Mean  optimal 
elasticity 
contrast  value 

0.1 

2.85  ±  0.0382 

5.62  ±  0.421 

0.2 

10.1  ±  0.367 

5.70  ±0.588 

0.5 

60.1 ±  4.19 

2.36  ±  0.393 

1.0 

80.2  ±  0.561 

2.47  ±  0.266 

2.0 

104  ±  3.42 

2.17  ±  0.422 
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Figure  3.  Three  candidate  automated  methods  for  MIE  boundary  condition  generation.  Top  row,  from  left  to  right:  surface  deformations 
calculated  from  diffusion  energy  matching,  Laplace  solution  energy,  and  thin-plate  spline  interpolation.  Middle  row:  target  registration  error 
(TRE)  distribution  for  each  method  when  compared  against  the  gold  standard  of  known  correspondence.  Bottom  row:  mappings  of  objective 
function  value  vs.  elasticity  contrast  ratio  (tumor:breast)  as  affected  by  the  boundary  condition  sets  generated  from  each  each  method.  The 
minimum  value  of  each  curve  corresponds  to  the  optimal  elasticity  contrast  that  can  be  achieved  by  the  algorithm  when  constrained  by  the 
inaccuracies  of  the  methods.  The  diffusion-based  mesh  is  both  qualitatively  and  quantitatively  the  worst  performer.  The  Laplace  solution 
appears  to  capture  the  shape  of  the  bladder  indentation  more  precisely,  but  the  thin-plate  spline  has  the  best  overall  accuracy  in  determining 
point  correspondence. 


Table  2.  Comparison  of  automated  point  correspondence  methods  on  MIE  reconstruction  quality.  Predicted  values  are  found  from  the 
minimum  point  of  the  curves  shown  in  Figure  3.  The  ideal  (true)  elasticity  contrast  is  6: 1 . 


Method 

Predicted  minimum 
objective  function 
value 

Reconstructed 
minimum  objective 
function  value 

Predicted 
optimal  elasticity 
contrast 

Reconstructed 
optimal  elasticity 
contrast 

Diffusion 

58.5 

58.8 

30.0 

12.6 

Laplace 

5.59 

5.59 

9.55 

10.3 

Thin-plate 

spline 

12.3 

12.3 

5.55 

5.66 

DISCUSSION 


The  results  of  the  boundary  condition  noise  experiment  indicate  that  improper  localization  of  boundary  points  greater  than  or 
equal  to  0.5  units  of  voxel  spacing  can  introduce  significant  error  to  the  reconstruction  process  and  impair  its  ability  to 
characterize  the  underlying  elasticity  distribution.  This  is  a  similar  result  to  prior  work  done  in  two-dimensional  systems  in 
which  successful  reconstructions  correlated  to  boundary  condition  selection  error  limited  to  half  a  pixel  length  [2].  It  also 
confirms  that  randomizing  the  vectors  is  a  significant  challenge  to  the  algorithm  because  it  introduces  highly  non-physical 
deformations  that  cause  backlash  in  the  finite  element  mesh  and  other  numerical  anomalies. 

Simulation  experiments  utilizing  automated  surface  registration  and  point  correspondence  methods  may  then  be  evaluated 
in  the  context  of  this  new  criterion.  Energy  matching  from  the  solutions  of  the  diffusion  and  Laplace  equations  yield 
boundary  condition  sets  that  are  inadequate  for  reconstructing  a  proper  elasticity  contrast.  This  can  be  partly  explained 
because  the  mean  errors  of  those  surface  registration  techniques  (as  compared  to  the  gold  standard  and  in  an  equivalent 
sense  to  the  tested  noise  levels)  are  approximately  3.3  and  1.7  voxel  units,  respectively,  and  are  far  too  large  for  the 
algorithm  to  handle.  The  diffusion-based  boundary  conditions  also  proved  more  difficult  to  obtain  a  stable  solution  for  in  the 
model,  which  probably  further  contributed  to  the  mismatch  in  reconstructed  elasticity  contrast.  However,  the  results  obtained 
using  the  thin-plate  spline  method  are  encouraging  because  the  mean  error  was  0.43  voxel  units,  thereby  satisfying  the 
threshold  while  demonstrating  reconstruction  success.  The  reconstruction  behavior  in  that  case  was  consistent  with  the 
predicted  objective  function  space  and  the  optimal  elasticity  contrast  was  found  to  be  within  6%  of  the  true  value.  This 
preliminary  result  appears  to  identify  the  use  of  thin-plate  spline  interpolation  as  a  strong  candidate  for  generating  boundary 
conditions  for  MIE.  The  implementation  for  that  method  required  the  use  of  40  control  points,  which  is  seen  as  a  reasonable 
choice  in  placing  fiducials  for  data  acquisition  in  order  to  capture  the  extent  of  anticipated  deformation  processes. 
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ABSTRACT 

Recent  advances  in  breast  cancer  imaging  have  generated  new  ways  to  characterize  the  disease.  Many  analysis 
techniques  require  a  method  for  determining  correspondence  between  a  pendant  breast  surface  before  and  after  a 
deformation.  In  this  paper,  an  automated  point  correspondence  method  that  uses  the  surface  Laplacian  or  the  diffusion 
equation  coupled  to  an  isocontour  matching  and  interpolation  scheme  are  presented.  This  method  is  compared  to  a  TPS 
interpolation  of  surface  displacements  tracked  by  fiducial  markers.  The  correspondence  methods  are  tested  on  two 
realistic  finite  element  simulations  of  a  breast  deformation  and  on  a  breast  phantom.  The  Laplace  correspondence 
method  resulted  in  a  mean  TRE  ranging  from  1.0  to  7.7  mm  for  deformations  ranging  from  13  to  33  mm,  outperforming 
the  diffusion  method.  The  TPS  method,  in  part  because  it  utilizes  fiducial  information,  performed  better  than  the 
Laplace  method,  with  mean  TRE  ranging  from  0.3  to  1.9  mm  for  the  same  range  of  deformations.  The  Laplace  and  TPS 
methods  have  the  potential  to  be  used  by  analyses  requiring  point  correspondence  between  deforming  surfaces. 

Keywords:  Registration,  non-rigid,  breast,  deformation,  correspondence,  surface  Laplacian,  diffusion,  finite  element, 
thin-plate  spline,  interpolation 


1.  INTRODUCTION 

As  breast  cancer  is  estimated  to  kill  over  40,000  women  and  be  diagnosed  in  more  than  178,000  in  2007  [1],  the 
detection  and  treatment  of  breast  cancer  is  an  important  area  of  scientific  research.  Many  novel  techniques  to  aid  in 
tumor  detection  are  being  developed  that  exploit  the  difference  in  physical  properties  between  healthy  and  cancerous 
tissue.  Some  of  these  techniques  measure  the  optical,  electrical,  or  elastic  properties  of  tissue,  including  near-infrared 
tomography  [6],  electrical  impendence  tomography  (EIT)  [7],  ultrasound  elastography  (USE)  [8],  magnetic  resonance 
elastography  (MRE)  [9],  and  in  particular,  modality-independent  elastography  (MIE)  [2,3]. 

MIE  is  a  reconstruction  algorithm  for  elasticity  imaging  that  can  be  used  for  detecting  breast  tumors.  It  involves 
imaging  a  pendent  breast  before  and  after  a  compression  and  using  these  images  to  reconstruct  the  elastic  properties  of 
the  tissue  using  a  nonlinear  optimization  framework,  computer  models  of  soft-tissue  deformation,  and  standard  measures 
of  image  similarity.  Unique  to  MIE  is  its  ability  to  utilize  images  from  any  modality  such  as  MRI  or  CT,  as  well  as  its 
usage  of  image  similarity  measures  that  make  direct  displacement  measurements  unnecessary. 

One  requirement  of  MIE  is  an  automated  method  of  finding  point  correspondence  between  the  pendent  breast  surfaces 
before  and  after  compression,  needed  to  specify  the  boundary  conditions  for  the  elasticity  model.  As  the  breast  is 
composed  of  soft  tissue  that  deforms  non-rigidly,  standard  rigid  registration  methods  cannot  be  applied.  Previous  work 
in  non-rigid  registration  includes  using  splines  and  FEM  models  [11],  as  well  as  point-based  methods  such  as  the 
symmetric  closest  point  (SCP)  algorithm  [10]. 

In  this  paper,  two  automated  methods  that  use  the  Laplace  and  diffusion  equations  to  establish  point  correspondence 
between  deformed  breast  surfaces  were  developed  and  compared  to  a  standard  thin-plate  spline  (TPS)  interpolation 
method  [4]. 
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2.  METHODOLOGY 

2.1  Laplace  and  diffusion  methods  of  finding  point  correspondence 

A  major  investigative  task  of  this  work  was  to  evaluate  whether  the  energy  distributions  modeled  by  a  partial  differential 
equation  (PDE)  over  an  undeformed  {source)  surface  and  a  deformed  {target)  surface  can  be  used  to  find  the 
correspondence  between  the  two  surfaces.  In  this  method,  the  Laplace  and  diffusion  equations  were  independently 
solved  over  the  source  and  target  meshes  using  the  finite  element  method  (FEM).  Laplace’s  equation  is  most  commonly 
used  to  describe  potential  flow  problems  such  as  in  thermal,  fluid,  and  electrostatic  systems  and  is  given  by 

■  q  (i) 

where  O  represents  the  potential  and  G  describes  the  spatially  varying  conductivity.  The  diffusion  equation  which 
allows  a  time-varying  potential  is  given  by 


g  =  (2) 

where  ®  represents  the  potential  and  a  is  the  diffusion  coefficient.  Let  ®S0Urce  refer  to  the  solution  to  the  Laplace  or 
diffusion  equation  over  the  source  surface,  and  let  ®target  refer  to  the  solution  over  the  target.  The  basic  premise  is  that 
the  potential  field  distributed  over  the  source  and  target  surfaces  as  calculated  by  the  Laplace  or  diffusion  equation  will 
provide  information  about  the  correspondence  between  the  source  and  target  surfaces. 

To  solve  the  equations,  Dirichlet  boundary  conditions  were  set  to  simulate  potential  flow  from  the  nipple  area  to  the 
chest  wall  over  the  surface  of  a  pendent  breast  (specifically,  nodes  in  the  nipple  and  chestwall  area  were  given  boundary 
values  of  1  and  0,  respectively).  To  solve  equation  1,  a  Galerkin  finite  element  method  is  used  whereby  the  equations 
are  expressed  along  the  surface  orientation  (a=l).  To  solve  equation  2,  a  similar  scheme  was  used  for  handling  the 
spatial  component  of  the  PDE  and  a  fully  implicit  backwards  Euler  scheme  was  used  for  time-stepping.  In  the  case  of 
equation  2,  a  no-flux  condition  was  prescribed  at  the  chest  wall,  and  the  potential  field  was  allowed  to  propagate  from 
the  nipple  (a=l).  In  this  calculation,  time-stepping  was  stopped  once  the  potential  field  reached  the  chest  wall. 

After  the  Laplace  or  diffusion  equation  was  solved  over  the  source  and  target  surfaces,  the  solutions  were  used  to 
establish  correspondence  between  the  source  and  target  nodes.  This  involved  two  distinct  processes:  finding 
correspondence  between  isocontours  of  ®source  and  ®target  and  then  “interpolating”  that  correspondence  back  to  every 
source  node  on  the  mesh.  In  the  first  step,  isocontours  were  extracted  from  ®source  and  ®target  for  a  set  of  selected 
isovalues.  The  correspondence  between  the  source  and  target  isocontours  was  determined  by  aligning  the  contours  by 
their  centroids  and  using  the  SCP  algorithm.  In  the  second  step,  the  displacement  vectors  at  the  source  isocontours  were 
interpolated  to  all  source  nodes  using  a  thin-plate  spline.  The  final  correspondence  was  found  by  adding  these 
displacements  to  the  source  nodes  to  get  the  location  of  the  corresponding  point  on  the  target  surface. 

The  method  can  be  summarized  in  the  following  steps: 

1.  Obtain  the  undeformed  source  mesh  and  deformed  target  mesh  that  define  a  breast  surface  before  and  after 
deformation. 

2.  Assign  boundary  conditions  at  nipple  and/or  chest  wall  nodes. 

3.  Solve  PDE  (diffusion  or  Laplace)  over  the  source  and  target  meshes  using  FEM. 

4.  Extract  isocontours  on  the  source  and  target  surfaces. 

5.  Determine  point  correspondence  between  source  and  target  isocontours  using  SCP. 

6.  Interpolate  displacements  at  source  isocontours  to  all  source  nodes. 

2.2  Using  thin-plate  spline  interpolation  to  find  point  correspondence 

One  advantage  of  the  PDE-based  correspondence  methods  is  that  they  do  not  explicitly  rely  on  external  markers  to 
constrain  the  matching  process.  However,  when  real-world  data  is  acquired,  fiducials  are  anticipated  to  be  available. 
Therefore,  TPS  interpolation  is  another  method  that  can  be  used  to  find  point  correspondence  [4].  Although  there  are 
many  different  methods  for  interpolation,  including  polynomial  splines  and  B-splines  [11],  TPS  interpolation  was  chosen 
in  part  because  it  does  not  require  a  regular  grid,  the  effects  of  changing  a  control  point  are  localized,  and  it  is  a  standard 


52 


method  that  has  been  successfully  used  in  many  non-rigid  registration  applications.  In  the  simulation  experiments 
described  below,  TPS  interpolation  was  used  to  find  point  correspondence  by  choosing  40  points  on  the  source  surface  to 
act  as  fiducials.  The  known  displacements  at  these  nodes  were  then  interpolated  to  all  surface  nodes  using  TPS. 

The  Laplace,  diffusion,  and  TPS  methods  for  finding  point  correspondence  described  above  were  tested  on  two 
simulation  data  sets  and  a  breast  phantom. 

2.3  Simulation  experiments 

To  perform  a  controlled  test  of  the  methods  described  above  on  a  breast-shaped  surface,  a  CT  image  volume  of  a 
pendant  breast  (courtesy  of  the  Dept,  of  Radiology,  University  of  Califomia-Davis)  was  segmented  to  create  a  source 
surface  consisting  of  6,313  points.  Two  types  of  deformations  were  simulated  by  assuming  different  contact  geometries 
of  an  air  bladder  being  inflated  against  the  breast  surface.  Circular  and  rectangular  cross-sectional  areas  of  a  Gaussian 
stress  distribution  positioned  along  the  lateral  aspect  of  the  breast  were  used  to  define  Type  2  boundary  conditions  for  a 
finite  element-based  deformation;  the  base  was  made  to  be  affixed  to  the  chest  wall.  The  displacement  solutions,  based 
on  a  three-dimensional  linear  elastic  model  of  a  Hookean  solid,  were  applied  to  create  the  target  surfaces  for  the  two 
simulations  (Figure  1). 


Figure  1.  Breast  surface  point  sets,  (a)  Source  surface  extracted  from  CT  volume  of  a  breast,  (b)  Target  surface 
generated  from  first  simulation  using  circular  cross-section  of  a  Gaussian  stress  distribution,  (c)  Target 
surface  generated  from  second  simulation  using  rectangular  cross-section  of  a  Gaussian  stress 
distribution.  The  correspondence  between  the  source  and  target  surfaces  was  determined  using  the 
Laplace,  diffusion,  and  TPS  methods. 

2.4  Phantom  experiments 

A  breast  phantom  was  constructed  to  test  the  point  correspondence  methods  with  real-world  data.  The  phantom  was 
fabricated  from  an  8%  w/v  solution  of  polyvinyl  alcohol  that  was  frozen  in  the  upper  half  of  a  2-liter  beverage  container 
for  16  hours.  After  8  hours  of  thawing,  thirty-four  1-mm  stainless  steel  ball  bearings  were  implanted  directly  under  the 
surface  of  the  resulting  cryogel  to  act  as  fiducials. 

The  phantom  was  then  imaged  inside  a  custom-built  rectangular  chamber  designed  to  deliver  compression  by  means  of 
an  air  bladder  placed  against  the  surface  of  the  phantom  (Figure  2).  CT  images  (512  x  512  x  174,  0.54  x  0.54  x  1  mm 
voxel  spacing)  were  acquired  with  the  phantom  at  three  different  states  of  mechanical  deformation  (undeformed,  50%  of 
maximum  bladder  pressure,  and  full  inflation).  Triangular  surface  meshes  were  obtained  by  semi-automatic 
segmentation  of  the  image  volumes  using  the  surface  extraction  tools  in  ANALYZE  6.0  (Mayo  Clinic,  Rochester,  MN), 
and  the  coordinates  of  the  fiducial  centroids  were  localized.  These  meshes  contained  approximately  8127,  6777,  and 
8260  nodes,  respectively. 
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Figure  2.  Experimental  system  for  image  data  acquisition.  A  polyvinyl  alcohol  cryogel  is  placed  within  a 
Plexiglas  chamber  with  its  surfaces  held  in  place  against  the  walls.  Compression  is  delivered  through  an 
air  bladder  (arrow)  inflated  manually  through  a  bulb  adapted  from  a  standard  sphygmomanometer. 

The  Laplace,  diffusion,  and  TPS  methods  were  tested  on  the  phantom  surface  meshes.  For  the  TPS  method,  30  of  the 
fiducials  was  used  in  the  interpolation  and  the  four  remaining  fiducials  were  reserved  for  validation.  The  fiducials  used 
in  interpolation  and  validation  were  selected  such  that  the  distribution  for  both  groups  over  the  surface  was  roughly  even 
and  included  the  deformed  region. 

2.5  Validation 

In  order  to  assess  the  accuracy  of  the  simulation  and  phantom  experiments,  the  target  registration  error  (TRE)  was 
calculated.  The  TRE  measures  the  error  between  the  correspondence  determined  by  the  registration  method  and  the  true 
correspondence  [11].  For  the  simulation  experiments,  the  TRE  was  calculated  as  the  Euclidean  distance  between  the 
corresponding  target  points  determined  by  the  Laplace,  diffusion,  or  TPS  method  and  the  true  target  points.  Since  the 
true  correspondence  between  the  source  and  target  surfaces  was  known,  the  TRE  was  calculated  for  each  source  node, 
and  the  average  and  maximum  were  reported.  For  the  phantom  experiment,  the  TRE  was  calculated  using  the  centroids 
of  the  bead  fiducials  implanted  directly  under  the  phantom  surface.  Since  the  “gold  standard”  correspondence  was 
known  only  at  the  fiducials,  the  TRE  could  only  be  calculated  at  these  locations. 

In  addition,  since  one  crucial  step  in  both  the  Laplace  and  diffusion  methods  is  to  find  point  correspondence  between  the 
source  and  target  isocontours  (step  6  of  algorithm  summary),  we  evaluated  how  well  the  SCP  algorithm  performed  in 
this  step  for  the  simulation  data.  To  accomplish  this,  the  SCP  method  was  given  a  set  of  source  isocontours  and  their 
true  corresponding  target  contours,  and  the  error  (the  Euclidean  distance  between  the  true  target  point  and  the 
corresponding  target  point  determined  by  SCP)  was  calculated. 

3.  RESULTS 

3.1  Simulation  1  (circular  deformation  source) 

The  Laplace  and  diffusion  equations  were  solved  over  the  surfaces  generated  from  simulation  1  (cranial-caudal 
deformation  source  with  maximum  displacement  of  33  mm)  to  find  point  correspondence  between  the  source  and  target 
breast  surfaces.  For  comparison,  TPS  interpolation  using  40  simulated  fiducials  was  also  used  to  find  point 
correspondence.  The  accuracy  of  each  method  was  assessed  by  calculating  the  TRE  at  each  surface  node  (Figure  3). 
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Figure  3.  TRE  displayed  for  breast  simulation  1  (circular  deformation  source)  when  (a)  Laplace  equation, 

(b)  diffusion  equation,  and  (c)  TPS  interpolation  were  used  to  find  point  correspondence.  The  TPS 
method  resulted  in  the  lowest  error  overall  (mean  TRE  0.4  mm),  followed  by  the  Laplace  method 
(mean  TRE  2.3  mm)  and  diffusion  method  (max  TRE  4.5  mm).  The  highest  TRE  is  found  in  the 
deformed  region  when  the  Laplace  method  is  used  and  in  the  base  when  the  diffusion  method  is 
used. 

The  results  (Table  1)  indicated  that  the  Laplace  method  performed  more  accurately  overall  than  the  diffusion  method; 
however,  the  area  with  the  highest  amount  of  error  differed.  When  the  Laplace  method  was  used,  the  deformed  region 
had  the  highest  error,  whereas  when  the  diffusion  method  was  used,  the  area  farthest  from  the  diffusion  source  had  the 
highest  error.  (In  this  case,  since  the  diffusion  source  was  located  in  the  nipple  area,  the  highest  error  occurred  in  the 
chest  wall  region.)  The  TPS  interpolation  had  the  lowest  error  overall,  and  the  error  distribution  over  the  surface  was 
related  to  the  locations  of  the  simulated  fiducials. 

The  results  given  above  pertain  to  a  simulated  compression  with  a  maximum  displacement  of  approximately  33  mm. 
Since  this  amount  of  compression  may  be  larger  than  is  needed  for  many  applications  and  may  introduce  other  unwanted 
effects  for  MIE  due  to  non-linear  elastic  behavior,  the  point  correspondence  methods  were  also  tested  for  lesser  amounts 
of  compression.  The  TRE  for  different  amounts  of  compression  when  the  Laplace  method  was  used  to  find  point 
correspondence  is  shown  in  Figure  4.  The  TRE  appears  in  Figure  2  and  increases  linearly  with  more  compression. 

The  mean  and  maximum  error  for  the  isocontour  point  correspondence  determined  by  the  SCP  algorithm  (detailed  in 
methods  section)  was  calculated  (Figure  5).  The  isocontour  correspondence  given  by  the  SCP  algorithm  had  a 
maximum  error  of  about  5  mm  for  the  maximum  compression  of  33  m. 
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Figure  2.  Maximum  TRE  (solid  line)  and  mean  TRE  (dotted  line)  when  the  Laplace  method  was  used  to  find  point 
correspondence  between  the  source  and  target  surfaces  at  different  levels  of  compression.  A  deformation 
of  100%  indicates  a  maximum  displacement  of  approximately  33  mm.  The  TRE  seems  to  increase 
linearly  with  respect  to  the  amount  of  deformation 


Fraction  of  deformation 

Figure  5.  Evaluation  of  the  accuracy  of  the  SCP  method  used  to  find  isocontour  correspondence.  The  mean  error  of 
the  SCP  method  (dotted  line)  is  compared  to  the  mean  error  of  the  Laplace  method  (solid  line)  when  SCP 
method  was  used  to  find  isoncontour  correspondence  for  simulation  1  data  at  different  levels  of 
compression.  (A  deformation  of  100%  indicates  a  maximum  displacement  of  approximately  33mm.) 
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3.2  Simulation  2  (rectangular  deformation  source) 

When  the  Laplace,  diffusion,  and  TPS  interpolation  methods  were  used  to  find  point  correspondence  between  the  breast 
surfaces  generated  by  simulation  2  (a  more  realistic  simulation  using  a  rectangular  deformation  source  with  a  maximum 
displacement  of  13  mm),  the  results  (Figure  6)  were  very  similar  to  those  from  simulation  1.  However,  the  TRE  for  all 
three  methods  (Table  1)  was  slightly  lower,  possibly  due  to  the  lower  degree  of  compression  simulated. 


meters 


(C) 

Figure  3.  TRE  displayed  for  breast  simulation  2  (rectangular  deformation  source)  when  (a)  Laplace  equation, 
(b)  diffusion  equation,  and  (c)  TPS  interpolation  were  used  to  find  point  correspondence.  The  TPS 
method  resulted  in  the  lowest  error  overall  (mean  TRE  0.3  mm),  followed  by  the  Laplace  method 
(mean  TRE  1 .0  mm)  and  diffusion  method  (mean  TRE  2.1  mm).  The  results  are  similar  to  that  of 
simulation  1 
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3.3  Phantom 


The  Laplace  and  diffusion  methods  were  used  to  determine  point  correspondence  between  the  noncompressed  and 
compressed  surfaces  of  a  breast  phantom.  The  results  were  validated  by  calculating  the  TRE  at  34  fiducials  located 
directly  below  the  surface  of  the  phantom.  For  comparison,  TPS  was  used  to  interpolate  the  displacements  of  30 
fiducials  to  all  surface  nodes,  and  the  TRE  was  calculated  using  the  4  remaining  fiducials. 

The  results  for  a  50  and  100%  compression  (with  a  maximum  displacements  of  about  20mm  and  36  mm,  respectively) 
are  shown  in  Table  2.  As  in  the  simulations,  the  Laplace  method  performed  better  overall  than  the  diffusion  method  and 
had  lower  TRE.  The  TRE  for  the  TPS  interpolation  was  lower  than  that  for  the  Laplace  and  diffusion  methods,  but 
varied  with  the  number  and  locations  of  fiducials  used  in  the  interpolation. 


Figure  4.  Breast  phantom  surfaces  (a)  before  compression,  (b)  at  50%  compression  with  maximum 
displacement  of  .020  m,  and  (c)  at  100%  compression  with  maximum  displacement  of  .036  m.)  Lines 
indicate  isocontours  at  different  values  of  k.  Black  nodes  at  the  nipple  and  base  indicate  the  nodes 
assigned  boundary  values. 
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Table  1.  TRE  for  different  point  correspondence  methods  tested  on  breast  surfaces  generated  from  simulation  1 
(point  deformation  source  with  max  displacement  of  33  mm)  and  simulation  2  (rectangular 
deformation  source  with  max  displacement  of  13  mm)  breast  surfaces.  The  Laplace  method 
outperformed  the  diffusion  method,  while  the  TPS  method  performed  best  of  all. 


Simulation  1 

Simulation  2 

Max  TRE  (mm) 

Mean  TRE  (mm) 

Max  TRE  (mm) 

Mean  TRE  (mm) 

Laplace 

14.6 

2.3 

4.6 

1.0 

Diffusion 

24.2 

4.5 

10.3 

2.1 

TPS  (40  fiducials) 

7.6 

0.4 

2.6 

0.3 

Table  2.  TRE  for  different  point  correspondence  methods  tested  on  breast  phantom  at  approximately  50%  and 
100%  compression,  with  max  displacements  of  20  and  36  mm,  respectively.  The  Laplace  method 
outperformed  the  diffusion  method,  while  the  TPS  method  performed  best  of  all. 


50%  Compression 

1 00%  Compression 

Max  TRE  (mm) 

Mean  TRE  (mm) 

Max  TRE  (mm) 

Mean  TRE  (mm) 

Laplace 

8.1 

3.5 

16.4 

7.7 

Diffusion 

11.9 

3.9 

19.4 

7.9 

TPS  * 

1.5 

1.1 

2.6 

1.9 

*  TPS  interpolation  using  30  fiducials;  4  fiducials  were  used  to  calculate  TRE. 


4.  DISCUSSION 

Of  the  three  registration  methods  evaluated,  the  TPS  method  consistently  outperformed  the  Laplace  and  diffusion 
methods  and  had  the  lowest  error  for  both  the  simulation  and  phantom  experiments.  However,  it  is  important  to  note 
that  a  comparison  of  the  PDE-based  methods  and  the  TPS  method  is  not  entirely  equal  since  the  TPS  method  relies  on 
fiducial  information  that  the  Laplace  and  diffusion  methods  do  not  require.  The  performance  of  the  TPS  method  is 
dependent  on  both  the  number  and  placement  of  these  fiducials.  These  results  indicate  that  30-  40  fiducials  with  an  even 
distribution  over  the  surface  should  be  sufficient  to  register  surfaces  (with  13-33  mm  deformations)  with  mean  TRE 
ranging  from  0.3  to  1.9  mm.  Although  further  studies  are  needed  to  determine  the  optimal  number  and  placement  of 
fiducials,  experience  suggests  that  increasing  the  number  of  fiducials  in  the  areas  with  greatest  deformation  increases 
registration  accuracy,  and  conversely,  lowering  the  number  of  fiducials  in  those  areas  causes  a  significant  decrease  in 
accuracy. 

The  results  indicate  that  the  Laplace  method  is  a  useable  surface  registration  method.  Although  the  Laplace  method  did 
not  perform  as  accurately  as  the  TPS  method,  it  has  the  advantage  of  not  requiring  fiducial  information.  However,  one 
of  the  challenges  of  the  Laplace  method  is  determining  the  regions  to  which  boundary  conditions  are  assigned.  Accurate 
selection  of  these  regions  is  important  because  the  implicit  correspondence  between  these  regions  is  used  by  the  Laplace 
equation  to  obtain  the  correspondence  for  the  rest  of  the  surface.  Lor  these  studies,  the  nipple  region  and  the  chest  wall 
boundary  regions  were  selected  manually,  further  studies  may  be  needed  to  find  a  method  to  automate  the  selection  of 
the  boundary  regions  and  to  evaluate  how  error  in  the  selection  of  these  regions  affects  the  final  registration  error. 

Although  the  diffusion  method  does  have  certain  advantages  over  the  Laplace  and  TPS  registration  methods,  several 
problems  prevent  it  from  becoming  a  viable  surface  registration  technique.  The  advantages  of  the  diffusion  method  are 
that  the  correspondence  near  to  the  diffusion  source  (in  this  case,  the  nipple)  is  relatively  accurate.  In  addition,  the 
diffusion  method  only  requires  boundary  conditions  to  be  set  in  one  region  (in  this  case,  the  nipple),  unlike  the  Laplace 
method,  which  requires  boundary  conditions  at  two  regions  (nipple  and  chest  wall  base),  and  the  TPS  method,  which 
requires  multiple  points  of  constraint  (at  34  fiducials). 
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However,  the  diffusion  method  does  not  appear  to  be  an  effective  surface  registration  method  for  the  following  reasons: 
the  results  indicate  a  substantial  amount  of  error,  the  registration  and  resulting  error  are  highly  dependent  on  the 
diffusion  parameters  chosen  (time  step  and  final  time  in  particular),  and  the  diffusion  parameters  must  be  manually 
adjusted  for  each  different  surface  mesh  since  there  is  no  automated  method  to  find  the  optimal  diffusion  parameters. 
Since  the  diffusion  described  by  the  PDE  is  by  definition  a  non  steady-state  process,  an  optimal  registration  requires  that 
the  diffusion  front  should  travel  over  the  entire  surface  between  the  nipple  and  base  and  stop  at  the  base  in  order  to 
assure  correspondence  for  as  much  as  the  surface  as  possible.  If  the  parameters  are  chosen  such  that  the  diffusion  front 
does  not  reach  the  base,  the  correspondence  for  the  regions  not  reached  by  the  diffusion  front  cannot  be  constrained  and 
must  be  interpolated  from  the  displacements  of  the  surrounding  regions.  Conversely,  if  the  diffusion  front  travels  for  too 
long  a  time,  the  solution  over  the  surface  approaches  saturation,  resulting  in  a  flat  gradient  and  lack  of  isocontours  from 
which  to  establish  correspondence.  Various  modifications  to  the  diffusion  method  employing  curvature  information  and 
using  different  diffusion  coefficients  were  tested,  but  none  were  successful.  Therefore,  the  sensitivity  of  the  diffusion 
method  to  parameters  and  substantial  amount  of  error  may  prevent  the  diffusion  method  from  being  a  viable  surface 
registration  method. 

The  TRE  measured  for  each  registration  technique  is  not  only  dependent  on  the  factors  described  above,  but  also  on  the 
amount  of  deformation  of  the  target  surface.  The  results  suggest  that  the  TRE  increases  roughly  linearly  with  the 
amount  of  deformation.  Using  the  simulation  and  phantom  data  presented  here,  one  may  be  able  to  estimate  the  range  of 
error  expected  when  one  of  the  described  methods  is  used  to  register  breast  surfaces  with  a  particular  amount  of 
compression.  Conversely,  the  maximum  amount  of  compression  that  will  yield  a  registration  within  a  given  error  bound 
can  be  roughly  estimated.  For  the  purposes  of  MIE,  realistic  compressions  will  be  in  the  range  of  1-2  cm. 

Another  factor  related  to  the  amount  of  compression  is  the  distribution  of  TRE  over  the  surface.  The  TRE  was  not 
evenly  distributed;  rather,  the  TRE  in  the  areas  of  greatest  deformation  tended  to  be  higher  than  the  TRE  elsewhere. 
Therefore,  the  mean  TRE  is  not  necessarily  the  best  measure  of  the  TRE  over  the  surface;  the  max  TRE  may  reflect  the 
error  in  the  deformed  regions  more  accurately. 

In  addition  to  the  evaluation  of  the  three  registration  methods,  the  performance  of  the  SCP  algorithm  was  evaluated  since 
the  matching  of  the  isocontours  extracted  from  the  source  and  target  surface  is  a  crucial  step  of  the  PDE-based 
registration  methods.  The  results  indicate  that  the  amount  of  error  the  SCP  algorithm  contributes  to  the  Laplace  and 
diffusion  methods  is  less  significant  small  when  compared  to  the  total  TRE  (Figure  2)  but  is  not  negligible. 

In  comparison  to  previous  studies,  the  Laplace  method  outperformed  the  modified  SCP  method  implemented  by  Schuler, 
et  al.  The  data  generated  by  first  simulation  described  in  this  paper  was  also  used  to  test  the  modified  SCP  method,  and 
whereas  the  Laplace  method  had  a  maximum  error  of  14.6  mm  for  a  deformation  of  33  mm,  the  modified  SCP  method 
had  a  maximum  error  of  27.8  mm  [5]. 

MIE  is  one  application  that  may  use  the  registration  methods  described  in  this  paper,  in  this  case  to  determine  boundary 
conditions  for  its  elasticity  model.  Preliminary  studies  indicate  that  TPS  is  the  most  viable  registration  method,  the  error 
of  which  is  within  the  bounds  required  for  a  successful  elasticity  reconstruction  (approximately  0.3  mm).  The  mean 
error  for  the  Laplace  registration  method  exceeds  MIE’s  error  bounds,  and  although  the  target  boundary  conditions 
produced  the  Laplace  method  resulted  in  a  viable  mesh,  the  resulting  elasticity  reconstruction  contained  a  considerable 
amount  of  error.  The  diffusion  method  could  not  be  used  in  conjunction  with  MIE  because  of  the  extreme  distortion  of 
the  target  finite  element  mesh  generated  from  the  surface  registration. 

5.  CONCLUSION 

The  results  of  the  simulation  and  phantom  experiments  indicate  that  while  TPS  interpolation  is  the  most  accurate  surface 
registration  method  of  those  evaluated,  the  Laplace  method  is  a  viable  surface  registration  technique  if  fiducials  are  not 
available.  Although  the  TPS  method  consistently  outperformed  the  Laplace  method,  its  performance  is  dependent  on  the 
number  and  distribution  of  fiducials  available.  Both  the  Laplace  and  TPS  methods  have  been  used  in  MIE  to  register 
breast  surfaces  in  order  to  determine  boundary  conditions  for  its  elasticity  model.  In  addition  to  MIE,  the  Laplace  and 
TPS  methods  also  have  potential  to  be  used  for  non-rigid  registration  in  more  general  applications. 
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Laplace/Diffusion  Registration  Method 

In  this  method,  the  potential  energy  field  modeled  by  either  the  Laplace  or  diffusion  equation  over  an  undeformed 
source  surface  and  a  deformed  target  surface  is  used  to  register  the  two  surfaces. 

1.  Obtain  undeformed  source  mesh  and  deformed  target  mesh  that  define  breast  surface  before  and  after 
deformation. 

2.  Solve  the  diffusion  or  Laplace  equation  over  the  source  and  target  meshes  using  FEM.  Laplace's  equation  is  used 
to  describe  potential  flow  in  thermal,  fluid,  and  electrostatic  systems  and  is  given  by: 

V  •  (oVO)  =  0 

The  diffusion  equation  allows  for  time-varying  potential  and  is  given  by 
d® 

—  =  V  •  (ctVO) 

where  CD  represents  the  potential,  o  the  conductivity,  and  a  the  diffusion  coefficient.  Use  following  type  1 
boundary  conditions:  set  0=1  at  nipple  nodes  for  both  the  Laplace  and  diffusion  equations,  and  set  O  =0  at 
chest  wall  for  the  Laplace  equation. 

3.  Extract  isocontours  on  the  source  and  target  surfaces. 

4.  Determine  point  correspondence  between  source  and  target  isocontours  using  Symmetric  Closest  Point  (SCP). 

5.  Interpolate  displacements  at  source  isocontours  to  all  source  nodes. 

6.  Assess  accuracy  by  calculating  the  target  registration  error  (TRE).  (For  the  simulations  described  below,  the  TRE 
was  calculated  at  each  surface  node.  For  the  phantom,  the  TRE  was  only  calculated  at  the  34  fiducials.) 


Obtain  source  and  target  mesh  and  set  boundary  conditions 


■  0  =  1 

■  0  =  0  (Laplace 

equation  only) 

source  surface 

target  surface 

Thin-Plate  Spline  (TPS)  Registration  Method 

TPS  interpolation  was  used  to  register  the  two  surfaces  by  interpolating  the  displacements  tracked  by  control  points 
to  all  surface  nodes.  The  displacements  at  each  surface  node  were  then  used  to  establish  the  correspondence 
between  the  two  surfaces.  For  the  simulations  described  below,  the  displacements  were  interpolated  using  40 
evenly  distributed  surface  nodes,  and  the  TRE  was  calculated  at  each  surface  node.  For  the  phantom 
experiment,  the  displacements  at  30  fiducials  were  interpolated,  and  the  TRE  was  calculated  at  the  remaining  4 
fiducials. 


EXPERIMENTAL  SETUP 


Simulation  Experiments 

To  test  these  methods,  a  CT  image  volume  of  a  pendant  breast  was  segmented  to  create  a  source  surface.  The 
breast  was  modeled  as  3D  linearly  elastic,  Hookean  solid,  and  two  types  of  deformations  were  simulated  by 
assuming  different  contact  geometries  of  an  air  bladder  being  inflated  against  breast  surface.  The  first  simulation 
used  a  circular  cross-section  of  a  Gaussian  stress  distribution  while  the  second  used  a  rectangular  cross-section;  the 
chest  wall  was  fixed  for  both.  The  resulting  target  surface  served  as  the  gold  standard  for  error  evaluation. 

Phantom  Experiments 

A  phantom  was  fabricated  from  a  polyvinyl  alcohol  cryogel  and  34  1- 
mm  stainless  steel  ball  bearings  were  implanted  directly  under  the 
surface  to  serve  as  fiducials.  The  phantom  was  imaged  by  CT  inside  a 
custom-built  rectangular  chamber  designed  to  deliver  compression 
through  an  air  bladder  (Figure  1),  at  three  different  levels  of 
deformation  (undeformed,  at  50%,  and  100%  max  inflation). 

Triangular  surface  meshes  and  coordinates  of  the  fiducial  centroids  Figure  l.  Compression  chamber  for  breast 
were  obtained  and  used  to  test  the  registration  methods.  phantom.  Arrow  indicates  inflation  bladder. 


Extract  isocontours  &  find  correspondence 
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Figure  2.  Accuracy  of  the  Laplace,  diffusion,  and  TPS  methods  when  used  to  register  breast  surfaces  from 
the  simulation  and  phantom  experiments  at  different  levels  of  compression.  The  TPS  method  consistently 
had  the  lowest  amount  of  error,  followed  by  the  Laplace,  and  the  diffusion  method  had  the  highest  amount 
of  error.  However,  note  that  the  TPS  method  cannot  be  directly  compared  to  the  Laplace  and  diffusion 
methods  because  it  utilized  fiducial  information  than  the  other  methods  do  not  require. 


Simulation  1  (circular  deformation  source,  max  displacement  33  mm) 


Figure  3.  Error  distributions  (measured  by  TRE)  when  the  Laplace,  diffusion,  and  TPS  methods  are  used  to 
register  breast  surfaces  from  simulation  1  (top  row)  and  simualtion  2  (bottom  row).  The  error  distributions  for 
simulation  1  and  2  are  similar.  For  the  Laplace  method,  the  highest  error  occurs  in  the  deformed  region;  for  the 
diffusion  method,  the  highest  error  occurs  near  the  base.  The  error  distribution  for  the  TPS  method  varies 
depending  on  the  placement  of  the  fiducials. 


Figure  4.  Effects  of  different  levels  of  compression  on  the  accuracy  of 
Laplace  registration.  The  max  and  mean  TRE  when  the  Laplace  method  is 
used  to  register  breast  surfaces  from  simulation  1  at  different  levels  of 
compression  are  plotted.  The  amount  of  TRE  seems  to  grow  linearly  with 
respect  to  the  amount  of  compression. 


CONCLUSION 


The  results  of  the  simulation  and  phantom  experiments  indicate  that  while  TPS  interpolation  is  the  most  accurate 
surface  registration  method  of  those  evaluated,  the  Laplace  method  may  be  a  viable  non-rigid  surface  registration 
technique  if  fiducials  are  not  available.  The  diffusion  method  does  not  seem  to  be  a  usable  registration  method  due 
to  the  high  error  and  the  difficulty  of  choosing  diffusion  parameters.  Although  the  TPS  method  consistently 
outperformed  the  Laplace  method,  its  performance  is  dependent  on  the  number  and  distribution  of  fiducials 
available.  Both  the  Laplace  and  TPS  methods  have  been  used  in  MIE  to  register  breast  surfaces  in  order  to 
determine  boundary  conditions  for  its  elasticity  model. 
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